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ABSTRACT 


The calculation of non-lifting potential flow about arbitrary three 
dimensional bodies is examined in detail with specific interest in the XYZ 
Potential Flow program developed by the David W. Taylor Naval Ship 
Research and Development Center. The program uses a surface singularity 
distribution to solve the Neumann boundary value problem by means of a 
source panel method assuming a flat element with a constant source 
density over the area of the element. Boundary conditions are applied at 
control points onthe elements producing a system of linear equations for 
the source density. When the source density is known, velocities and 


pressure coefficients may be calculated. 


The main purpose of this paper is to present the details of the 
approximation of an arbitrary three dimensional body using quadrilateral 
elements, and to provide a detailed derivation of the exact source panel 
Integrations in order to gain insight for future research at Texas A&M 
University. A variation of the Hess method of surface discretization 
using quadrilateral source panels is described in detail as it is used in the 
XYZ Potential Flow program. The exact source panel integrations are 


derived in detail. 


A general discussion of other aspects of the program is included. 
Velocities and pressure coefficients for flow about a triaxial ellipsoid 
are calculated using the XYZ Potential Flow Program solution, and the 
results are comparedwith the analytical solution and the Hess Program 


solution. 
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1.0 INTRODUCTION 


This paper examines two aspects of the development of the xYZ 
Potential Flow Program (hereafter referred to as the XYZPF Program), a 
FORTRAN program which uses a source panel method to approximate 
solutions to steady potential flow problems about arbitrary three 
dimensional bodies. The aspects examined in detail are (1) the 
description of the details of the approximation of an arbitrary 
three-dimensional body using quadrilateral elements, and (2) a detailed 


derivation of the exact source panel integrations. 


The XYZPF Program was developed specifically for applications in 
Numerical ship modelling and hydrodynamics studies at the David W. 
Taylor Naval Ship Research and Development Center (NSRDC) in Bethesda, 
Maryland. The format of the program is based onthe work of Hess and 
Smith (1962) in the numerical calculation of non-lifting potential flow. 
Asimilar program is maintained by the Aerodynamics Division of the 
McDonnell-Douglas Corporation, referred to in this paper as the "Hess 
program.” The XYZPF Program is a modification of what has come to be 
known generally as the Hess Method. The most significant modifications 
are improvements to the method of solving the matrix equation for the 
source density, and greater flexibility in the input options (Dawson and 
Dean 1972). 


Though potential flow is a product of mathematics, and is never 
found in a real fluid, the results of potential flow calculations provide 


usable information for flow regions external to a thin boundary layer, 





with little or no boundary layer separation. For such flow fields, the 
region outside the boundary layer may be considered to be effectively 
inviscid, and may be closely approximated by potential flow models. 
Small viscous effects canbe accountedfor by “thickening” the body by the 
appropriate displacement thickness. Displacement thickness accounts for 
the region of retarded fluid flow in the boundary layer inversely 
proportional to the square of the free stream velocity. Downstream of 
the point of boundary layer separation, the potential flow model no 


longer applies. 


Prior to the development of numerical methods, analytical 
solutions were generally restricted to simple analytical shapes (Kellogg 
1929). The need to solve boundary value problems for arbitrary 
boundaries in continuum mechanics has fostered the development of 
numerical approximations to the integral equation expressions. While 
the integration methods have been well known for quite some time, only 
since the advent of high speed computers have many of the problems 
been practical to solve by numerical methods. Amongthe numerical 
methods being used are finite differences, finite elements, and the 


boundary element method. 


"Finite differences” and "finite elements” are numerical methods 
which satisfy the boundary conditions, and then approximate the 
solution to the governing equation in the fluid domain. These methods 


discretize the domain into a network of elements or cells. 





Another approach is what is now known as the "Boundary Element 
Method,” in which the governing equation is exactly satisfied in the 
domain, and the boundary conditions are applied through a boundary 
discretization method. The boundary value problem is reformulated as a 
boundary integral equation which is then discretized by subdividing the 
boundary into a finite number of surface elements. Each element is 
represented by an analytical function, and the source density function is 
integrated over the surface of each element. Two factors governing the 
accuracyof the boundary element method are the boundary discretization 
method and the source panel integration. These two factors are examined 
in detail in this report, as a detailed derivation of the exact source panel 
integration, including the development of the source panel geometry, has 


not previously appeared in literature. 


The difference between the domain methods and the boundary 
methods is significant. The domain methods discretize the domain, while 
the boundary methods discretize the boundary. Thus, the boundary method 
reduces the dimension of the problem by one, as depicted in figure (1). In 
the application of the XYZPF Program, the problem is reduced from a 
three-dimensional problem in the domain to a two-dimensional problem 
on the boundaries. This method is well suited to problems in which the 
limits of the domain are infinite or difficult to define, in that the 


problem is applied to the boundary rather than the domain. 








FINITE ELEMENT DISCRETIZATION 


© 


BOUNDARY ELEMENT DISCRETIZATION 


Figure 1. Discretization Methods 


Just as there are many variations of domain methods, there are 
also a variety of boundary methods. In general, they canbe classified as 
“indirect” or “direct” formulations. The “indirect” method assumes a 
continuous source distribution over the surface of the body, and a 
solution which satisfies both the governing equation in the domain, and 
the boundary conditions onthe body surface. The result is an integral 
equation onthe boundary which has the surface source density function as 
its unknown. By enforcing the boundary conditions at control points on 
the surface, a system of equations is produced by which the source 


density may be determined. 


The “direct” method solves the velocity potential function through 
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an application of Green’s Second Identity requiring the solution of a 
source distribution and a dipole distribution onthe boundary. The direct 
method has more physical significance to the general boundary value 
problem, and more versatility in its application as it canbe applied to 
Neumann problems, Dirichlet problems, or mixed boundary value problems 


(Brebbia 1984). 


The simplicity and accuracy of the indirect method has made it 
attractive for many applications. The source panel method is an 
application of the indirect formulation of the boundary element method 
to the Neumann type of potential flow problem, for which the normal 


derivative of the potential function is prescribed onthe boundary. 


1.1 OBJECTIVES 


The purpose of this paper is (1) to describe the details of the 
approximation of an arbitrary three-dimensional body using quadrilateral 
elements, and (2) to provide a detailed derivation of the exact source 
panel integrations for use in future investigations at Texas A&M 
University using panels of higher order geometries and source density 
functions. This paper is not intended to be a user's manual, though a 
general discussion of other aspects of the program is also included. 
NSRDC Report 3892 (Dawson and Dean 1972) is a summary of the XYZPF 


Program for those strictly interested in its use. 





2.0 HISTORICAL DEVELOPMENT 


The foundations of the boundary element method were laid early in 
this century beginning with Fredholm in 1903 when he established the 
existence of solutions to the Neumann problem through a reconstruction 
of the problem using a discretized boundary (Kellogg 1929). The solution 
was determined to be the potential of a simple source distribution ona 
boundary with a continuous normal derivative for an infinite domain. 
Later works by Kellogg (1929) in potential theory demonstrated the 
application of the boundary integral equation method in electrostatics, 
heat transfer, flow in porous media, and fluid flow problems, but 
development was limited by the difficulty of obtaining analytical 
solutions. No significant advances were made until interest in boundary 
integral equation methods was revitalized with the advent of high speed 
electronic computers. Investigators were then able to discretize the 
boundaries and solve the integral equations numerically. This method of 
solution became known as the boundary element method. Early 
development of such numerical methods was pioneered by Hess and Smith 
(1962) and Jaswon and Symm (1963). Hess and Smith dealt primarily with 
the indirect formulation eventually leading to a solution for the three 
dimensional problem as described in this paper. In a parallel work, 
Jaswon and Symm developed a direct f ormulation approachto the two 
dimensional problem. The XYZPF Program is based primarily onthe work 
of Hess and Smith. Hess has since developed a higher order panel method 
(Hess 1979) and Lefebvre modified the XYZPF Program for calculating 


velocity potentials for five degrees of freedom (Lefebvre 1982). 





3.0 THEORETICAL DEVELOPMENT 


3.1 THE POTENTIAL FLOW PROBLEM IN THREE DIMENSIONS 


The governing equation for ideal (incompressible, inviscid, 


irrotational) flow is Laplace's equation: 


V29 =0 (1) 


where $ is the velocity potential, and V* is the Laplacian operator. The 
XYZPF Program deals with steady, uniform flow of an ideal fluid about an 
arbitrary three dimensional body. The velocity components at any point 
within the flow field may be obtained from the negative gradient of the 


velocity potential, i. e. 





The freestream flow V is defined as a uniform stream of unit 


magnitude. 


- 1 (2) 





The key to the boundary element method is the Divergence Theorem 
(Green’s Theorem) which relates a volume integral to an equivalent 


surface integral reducing the three-dimensional problem to a 
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two-dimensional one. The expression for Green's second identity is 


(Lamb 1924): 
I (Фууу - муф) 09 = ІІ War $ an ) dT (4) 


in which Q represents the integration over the three dimensional domain, 
and T represents integration over the two dimensional boundary. The 
partial derivatives are taken with respect to the outward normal, n. The 
weighting function, w, is usually chosento be the fundamental solution 
for three dimensions, w = 1/(4trr), where r is the distance from the 


source to an arbitrary point onthe boundary. 





Figure 2. Potential Flow in Three Dimensions 


Consider an arbitrary three-dimensional body with surface I, 
having an equation of the form F(x, y, z)= 0 where x, y, z are Cartesian 
coordinates of the global reference system as shown in Figure (2). The 
unit outward normal, n, at any point onthe surface is given by the 


gradient of the function describing the surface divided by the magnitude 
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of the gradient, i.e. 
+VF 


vr 





(3) 


where the sign of the unit normal vector is chosento ensure that the 
vector is an outward normal. The potential function $ describing the 


flow field must meet the following boundary conditions: 
a. V$ - 0 (Laplace's Equation) (6) 


b. For an impermeable boundary, the velocity normal to the surface 
must be zerorelative to the boundary (the Neumann boundary 


condition): 


aby 


* Qn *T 


c. The velocity potential approaches the freestream velocity 


potential as the distance from the body goes to infinity: 


e 55 г | Pg (8)‏ + چ 


The total potential at any point in the domain is composedof the 


freestream potential and the disturbance potential due to the body, 


Ses 8 (9) 





The disturbance potential, $, satisfies the following boundary 


conditions: 
а. У29 - 0 (10) 


b. From equation (7), the velocity normal to the boundary due to 
the disturbance and due to the onset flow must be of equal 
magnitude, but opposite sign. Then from equation (9) 

rans 


Note that the normal vector is a function of position onthe surface 


of the body. 


c. The disturbance potential approaches zero as the distance from 
the body goes to infinity, i. e. 


n 0) a5 Ir | — co (12) 


` 


3.2 MATHEMATICAL MODEL 


The disturbance potential of the body may be represented by a 
distribution of a source density function c over the body surface. The 
potential at an arbitrary point P(x, u, z) due to the surface potential is 


(Kellogg 1929): 





E eso (13) 


where q is the integration point onthe surface, and 





is the distance from the field point P to the integration point q. 


The source density distribution function must satisfy the boundary 
conditions for the disturbance potential. Boundary conditions (10) and 
(12) are automatically satisfied by the form of the integrand. However, 
equation (11), the velocity normal to the boundary, combined with the 
Neumann boundary condition, equation (7), is the key to solving the 


boundary integral problem. 


The integrand becomes singular as the surface of the body is 
approached, i. e. ir] goes to zero. The singularity represents the local 
fluid flux norma! to the boundary due to the local source density. The 
principal value of the singularity is -21ro(p), determined through a 
limiting process of the Gauss Flux Theorem (Kellogg 1929). The point p 
represents a field point which lies onthe boundary. The integral 
expression is now composedof the contribution of the local source 
density and the contribution of the source density function over the 
remainder of the body surface. Solving for the velocity normal to the 


surface yields the following expression: 





D a 14) 
CAR Rep E En Ear 14 


From equation (11), this expression becomes: 


zro{n) — ШЕ an E | а sun pv. Cim 


This equation is a two dimensional Fredholm integral equation of the 
secondkind, which ensures a unique solution, and that the diagonal 
elements of the system matrix will be dominant, each having a value of 
27 (Kellogg 1929). Once equation (15) has been solved for the source 
density o, the velocity components at any point of the flow field may be 
obtained by differentiating the disturbance potential function (13) with 


respect to the coordinate directions and adding the components of the 


freestream flow, Vo. 


t ач e M GN 
Міх, Ц, 2) = М – e - non) - x (10) 
c 


The body shape does not have to be slender, axisymmetric, or 
simply connected, allowing for analysis of interior flow and a wide range 
Of applications of the method. The only restriction imposed onthe form 
of the body is that it must have a continuous normal vector. 
Discontinuities in the right hand side of equation (15) will produce 
unwanted singularities. Thus, in the processof approximating a body 
which has distinct corners, where there is clearly a discontinuity in the 


normal vector, the corner must be replaced by a surface with some finite 
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where po, is the static pressure at infinity. 


3.3 NUMERICAL MODEL 


In order to represent the surface of a body in the domain 
mathematically, the body may be described by analytical expressions 
which may provide an exact representation of the surface. However, the 
types of bodies which canbe adequately described by such methods are 
severely limited. Another way to represent the body is to use a large 
number of analytical expressions, each describing only a small portion of 
the body. Hess and Smith (1962) suggested the use of an assembly of flat 
quadrilateral elements to model the actual surface of the body, as shown 
in Figure (3). Each quadrilateral approximates a region of the surface 
described by points which lie onthe actual surface of the body. As planar 
elements, these quadrilaterals are clearly analytical, and when carefully 
constructed, the elements can approximate arbitrary three dimensional 


body surfaces without restriction. 





Figure 3. Approximation of the body by surface elements 
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curvature. However, application of this method has shown that the flow 
calculations give correctresults for convex corners, while unrounded 
concavecorners may or may not produce significant error, depending on 


the angle produced by the corner (Hess and Smith 1962). 


Because of the method of approximation, the calculation of flow 
velocities onthe body surface are restricted to the points at which the 
boundary conditions were applied. Velocities at points other than those 
must be obtained by interpolation. Direct calculation of velocities at the 


edge of an element yields infinite velocities. 


With the solution of the system of linear equations for the source 
densities, the flow velocities at any point in the domain may be obtained 
from equation (16), and pressure coefficients are then computed from the 


velocities using a form of the Bernoulli equation: 


2 - ad y = 
ar ( Ir d 
EE 


Pit! = ЕЕ + 
р 





4 


NE 
where P(t) is a constant independent of position. In the XYZPF Program, 
the flow is steady. Therefore, equation (17) canbe reduced to 

pe > P V É = constant (18) 


and the pressure field canbe expressed in terms of a dimensionless 








pressure coefficient Co as: 
_ y | 
Eas د‎ т 
P 1 B 2 
72 P| Voc БА 
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The XYZPF Program uses the discretization procedure described by 
Hess and Smith (1962) with some minor modifications. The three 
dimensional body surface may be described using a large number of plane 
quadrilateral elements, each assumed to have a constant source density 
over the area of each element. Regions of the body requiring higher 
resolution for sharp curvature or anticipated velocity gradients will 


require a higher concentration of elements. 


Because the plane quadrilateral elements cannot fit edge to edge on 
a rounded surface, small gaps in the panel approximation contribute to the 
error of the approximation. However, the error due to the gaps is 
negligible when compared with the error of the basic model, that is, using 
flat panels to approximate a curved surface. Triangular elements have 
been suggested in an attempt to eliminate the gaps (Levu 1959), but the 
increased accuracyis so small that it may not justify the additional work 
of organizing the triangular elements in lieu of the simpler quadrilaterals 
(Hess and Smith 1966). The method presented is valid for an polygonal 


element with any number of sides. 


Equation (15) cannow be decomposedinto a summation of 
integrals, each representing the contribution of one element of the body 
surface. The unknown source density canbe taken outside the integral, 
since it is assumed to be a constant over each element. The integration is 
performed over the area of the source element, and the boundary 
condition equation (11) is then enforcedat a single point p in each 


remaining element. By performing this operation at each element of the 
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surface, a system of linear equations is generated which is equal in 
number to the number of surface elements and the number of unknown 
source densities. Equation (15) cannow be approximated by the matrix 


equation (Dawson and Dean 1972): 


Я — 20) 
O, - > Vj Cis + Y; ( 
1 
Where 
poa = [| on .. ПА 
y en an, Ar 
I 


| 
Ше =. Ay: M 


ha 


It is important to note that the influence coefficients Ci; and Ci; are 


functions of geometry only, and once computed, need not be recomputed 
for analysis of several different flows. From the solution of equation 
(20) onthe discretized surface, equation (13) may be applied at any point 
in the domain. Then, the velocity at an arbitrary field point P(x, u, z) in 
the domain may be determined from equation (16). With the velocity 


known, the pressure coefficient is determined from equation (19). 





4.0 ORGANIZATIONOF THE PROGRAM 


The XYZPF Program is actually composedof seven independent 
programs, referred to as sections PFI through PF7, each of which builds 
on data generated from a previous section. This type of organization 
allows the user the flexibility of rerunning portions of the program using 
different flow parameters without having to go through the time 
consuming process of recalculating the influence coefficient matrix, 
which is dependent only onthe geometry of the body. While the NSRDC 
program is very similar to the Hess program, there are also some 
significant differences. The following list of differences is taken from 
NSRDC Report 3892 (Dawson and Dean 1972): 


(1) The input to XYZ PF is arranged to facilitate the preparation of 
input for a series of problems in which only one part of the body is 
changed. Also, a number of checks are made onthe input to help detect 
errors. 

(2) Anoption was added for the recomputation of the source 
density and velocities for only part of the body when only small changes 
are made. This option also provides for the use of the solution of one 
problem as an initial guess for the solution of another problem. 

(3) The matrix of influence coefficients is computed column by 
column instead of row by row. This column arrangement was used for the 
original LARC computer version because it required much less high speed 
memory. The computation is also about 10% faster this way than with 
the row-by-row arrangement. 


(4) Asimultaneous displacement iterative scheme is used to solve 
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the matrix coefficient for the source density. The scheme is slower than 
the successive displacement (Gauss-Seidel) scheme used in the Hess 
program, but it canbe carried out using the matrix column by column 
instead of row by row. 

(5) When possible, an extrapolation procedure is used to reduce the 
number of iterations required for convergence. One such method is the 


Richardson extrapolation. 


The methods used in the XYZPF Program will be discussed in detail 


in the following sections. 





9.0 DETAILS OF THE SURFACE APPROXIMATION 
9.1 PREPARATIONOF THE INPUT 


Section PF1 is set up to read and store the input data, and to 
examine the cornerpoints of the quadrilaterals to detect obvious errors in 
the input. Because of the number of points which may have to be entered 
for a complex geometry, the user input is a major source of program 
error, and this first look for input errors will save lot of run time in the 
program as a whole. If Section PF! detects major errors in the input, the 
program will not continue with the calculation of the coefficient matrix, 
but will stop and identify the grid location of the error. Minor errors may 


not cause the program to stop, but will be noted in the output. 


One of the major advantages to this program is in the organization 
of the input data. The surface is input in sections so that small portions 
of the input geometry may be changed without having to recalculated the 
points for the entire surface. The program also takes advantage of 
symmetry to minimize the input effort. Only the portion of the body 
which has no redundancy needs to be entered point by point. The 
remainder of the body is reflected across the planes of symmetry by the 


program to complete the surface representation. 


The surface is represented by a set of points in three-dimensional 
Space which lie onthe actual surface, and which will later be used to 
define the plane quadrilateral sourceelements. These points are defined 


in the global reference system. The points on the surface should be 
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selected in such a way as to provide an accurate representation of the 
surface with the fewest number of points possible. Portions of the 
surface which are highly curved will require a larger number of points to 
provide adequate resolution. Additionally, portions of the surface in 
which the flow field is expected to change rapidly will require a large 
number of points to accurately determine the flow field in that region. 
Some familiarity with fluid dynamics will provide a somewhat intuitive 
approachto properly distributing the elements. Elements should change 
gradually in size from areas of high concentration to those of just a few 
large elements, changing no more than 50 percent in size between 
adjacent elements (Hess and Smith 1966). The accuracyis only as goodas 
that provided by the largest element in a particular area. The use of 


quadrilateral elements facilitates the use of known analytical equations 


and body contours to determine the input points. 


ORIGIN 
Figure 4. The 3D quadrilateral element in global coordinates 
For the purposes of this program, the body surface is represented 
by a large number of plane quadrilateral elements as shown in figure (3), 


each of which is assumed to have a constant value of source density over 


the area of the element. Eachelement is defined by four input points 


20 





which lie onthe actual surface as shown in fiqure (4). Since each input 
point can be used as a corner for up to four elements, there is no need to 
enter the same point four separate times. The input points are organized 
in groups of four to form the quadrilateral element, and each point may 
also be associated with adjacent quadrilaterals. This is accomplished 
through the use of a two dimensional coordinate system in which the user 
assigns a pair of integers, m and n, to each point which identifies the 
"row" and the "column" in which it lies. A column of points will be given 
a common value of n, and each point in that column will have a unique 
value of m corresponding to the row in which the point lies. The 
orientation of these "coordinate" integers determines the direction of the 
outward normal for each element. Looking from the flow field toward 
the section of elements, if the values of m are increasing upward, the 
values of n must increase to the right. Increasing m and n can point in any 
direction with respect to the global reference system. In fact, the 
orientation can change from one section to another. However, any other 
relationship between m and n will produce an incorrect normal vector. 
Once assigned, the values of m and n also serve to identify the element 
for which the correspondingpoint is a corner. The four points which 
form a quadrilateral element are two points of one column, or n line, with 
consecutive m numbers, and two points of the next higher n line with the 
same m numbers as the previous two points. Thus, the element m, nis 
composedof the points identified by (m, n), (m*1, n), (m, n*1), and (m-1, 


п+1). 
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Figure 3. Quadrilateral index numbers 


Each section of the body surface is formed by specifying a set of 
corner points corresponding to the m, n pairs for all of the quadrilaterals 
of the section. The user will sequentially assign an m number to the 
points for eachn line, and also number the n lines for the section points 
entered. The first point In each n line will always have m = 1. Then lines 
are also numbered sequentially, but the value of n is not reset for each 
new section. The sequence of n numbers runs through all the sections as 
shown in figure (5). Points ona particular row or column do not have to 
be strictly colinear. By forming nearly triangular elements, a rounded 
planform canbe approximated without conflicting with the numbering 


convention, as shown in figure (6). 
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Figure 6. Approximating a thin region with rounded planform (Hess 4 Smith 1962) 


By entering data in sections, small changes in geometry can be 
performed without having to reenter all the points associated with the 
body surface. This feature is unique to the XYZPF Program, and offers a 
great deal of flexibility in design work. However, with the added 
flexibility comes more restrictions on both the input of the original 
geometry and onany modifications. There are four important 
restrictions onthe input which are required to provide quadrilateral 
elements in groups of four to facilitate geometry calculations (Dawson 
and Dean 1972): 


(1) There must be an even number of elements in both the m andn 


directions in each section of the body. 


(2) The common corner point of a group of four elements must not 
coincide with any other corner point. The sides between the elements 
serve to define the local coordinate system, and serve as the axis of 


rotation when the surface is flattened for numerical differentiation of 
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the velocity potential. 


(3) Each set of four elements must have at least seven distinct 
corner points to allow the elements to more closely conformto a curved 
surface. This also allows for convergenceof N-lines or M-lines as may 
occur, for example, at the leading edge of an ellipsoid. Thus, only two of 
the four quadrilateral may degenerate into triangles by having two of 
their corner points coincide. This does not necessarily eliminate the 
possibility of more than two "triangular" elements since the adjacent 


sides of a quadrilateral may be colinear as shown in figure (6). 


(4) The normal vectors between two adjacent quadrilaterals ina 
group of four must be less than 90 degrees and preferably less than 45 
degrees. If a sharp edge is required, it should be a concavecorner with 
respect to the flow field, and the input should be arranged so that the 
edge is along an outside boundary of the groups of f our, and not through 


the center. 


When making small changes to the original geometry, the number of 
elements used in anew section must be the same as the number used in 
the original section unless the part being changed is at the end of the 
input data. Section configurations may be selected by natural divisions, 
as amatter of convenienceto more easily handle large numbers of points, 


or as a tool to take advantage of symmetry. 


In setting up input data to use planes of symmetry, it is important 


to note that the XYZPF Program has certain restrictions onthe choice of 
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symmetry planes. The user only has the option to select the number of 
symmetry planes. The planes which will be used as symmetry planes are 
preselected by the program to optimize the calculation procedure. 
Therefore, knowing this, the preparation of input data must consider the 


following restrictions imposed by the program (Dean and Dawson 1972). 


If only one plane of symmetry is used, the plane of symmetry is the 
y = O plane of the global coordinate system. Assuch, all the y 
coordinates of the input points must be of the same sign, i. e., all 
positive or all negative. If the body is closedand intersects the plane of 
symmetry, the points touching the plane, i. e., corresponding to у = 0, 


must also be entered with the input points. 


If two planes of symmetry are used, the planes of symmetry are the 
u = O plane and the z= O plane in the global coordinate system. Again, the 
y coordinates of all input points must have the same sign, positive or 
negative, and the z coordinates of all points must be of the same sign, 
positive or negative without regard to the sign of y. If the body surface 
intersects one or both of the planes of symmetry, the points which lie in 
the plane, i. e., those correspondingto y = 0 or z= 0, must also be entered 


with the input points. 


If three planes of symmetry are used, clearly the planes are the 
reference planes of the global coordinate system. As with the previous 
cases, all the x coordinates of the input points must be of the same sign, 
and similarly for the y and z coordinates. If any part of the body 


intersects any of the planes of symmetry, the points which lie in that 
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plane, i. e, x 0, yu- Oorz 7 0, must also be entered with the input 


points. 


9.2 SOURCE PANEL GEOMETRY 





Figure 7. The outer normal to the quadrilateral element 


With the surface points identified bu the location numbers, m and n, 
and arranged in accordancewith program requirements, calculation of 
various aspects of the source panel geometry and formation of the plane 
quadrilateral element is the next step in the numerical integration 
process. Formation of all of the planar elements is identical, so the 
following discussion of source panel geometry will deal with only one 
characteristic element. The four corner points forming the basic 
quadrilateral are numbered in a clockwise direction from I to 4 as shown 
in figure (7). It does not matter which corner point is identified with the 
number 1 subscript, but the remaining points must be numbered 
consecutively in a clockwise direction when observed from the flow field 
in order to ensure an outward directed normal vector. These subscripts 
will be used to identify the cornerpoints for the remainder of this 


discussion. For this example, the points will be identified as follows: 
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Position Numbers Global Coordinates 


m, n ^i, Yi, £ 

m*l n Ao , Yo , £o 
m+l, n+l AE 
m, n*1 Хд, Үд, Za 


In forming the plane quadrilateral elements, the corner points, 
which are generally not coplanar, are used to form the local coordinate 
system, relative to the element. Recalling that the crossproduct of two 
vectorsyields a vector solution which is perpendicular to both of the 
original vectors, the normal to the element may be obtained from the 


crossproduct of the diagonals of the element, 
N = Т4 x Тіз (21) 
where T;¡z is the vector from point 1 to point 3, and T54 is the vector 
from corner point 2 to point 4. The unit normal ts then: 
ls 


Ta a x hz "E 


This unit normal now represents the first of the three local 
coordinate directions, this one in the € direction. The side of the 
quadrilateral from point 2 to point 3 is then used to obtain the second 


coordinate vector. 
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Figure 8. The second local coordinate vector 
N, = n x LEX: (23) 


and the unit vector 


“% 
а= 


is = — 
SEIS 


(24) 


Similarly, the third local coordinate vector is obtained from the 


cross product of n» and n, the result of which is a unit vector. 


Ша 


noo 


Ne, X 
Figure 9. The third local coordinate vector 
ns = пх n (29) 


The unit vectors ns, n», and n form an orthonormal basis and define 


the local coordinate system for the element in the & 7, апа č directions 
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respectively. Other methods of obtaining an orthonormal basis could be 
used just as well, and would make no difference to the remaining 
computations. The origin of the local coordinate system would most 
correctlybe located at the "null point,” the point at which the velocity 
potential has no contribution to the tangential velocity component on the 
source element. The null point 15 the point in each quadrilateral element 
where the normal velocity boundary condition is applied. However, with 
the exception of long, thin quadrilaterals, the physical difference between 
the null point and the centroid of the quadrilateral is not significant. The 
AYZPF Program will print a warning in the output when a quadrilateral is 
long and thin enough to jeopardize the accuracy of the approximation in 
that region. By locating the origin of the local coordinate system at the 
centroid, rather than at the null point, the difficult processof locating 
the null point for each element canbe eliminated, later calculations of 
the multipole expansion canbe simplified, and the boundary conditions 
can be applied at the centroid without contributing significant error to 
the approximation (Hess and Smith 1966). Therefore, the origin for each 
local coordinate system is located at the centroid for the respective 


element. 
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9.3 LOCATING THE CENTROID 


| 


Figure 10. Locating the centroid of the quadrilateral 


The centroid of the element may be calculated by first dividing the 
area of the quadrilateral into two triangular areas, the triangles being 
separated by the line from point 2 to point 4. The area A, of the triangle 
defined by corner points 2, 3, and 4 is 


A = Tax Tos 


1 





(26) 
Similarly, the area A, of the triangle defined by corner points |, 2, and 4 
IS 


bo = y [ex Ta| (27) 


In the global coordinate system, the X component of the centroid is 


given bu 


- 


xc | 


= Aix; + Az a2 (28) 
A, + 57 
where Xand X> are the averages of the X components of the corner 


points of each triangle. Substituting the values for X, and Xo: 
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1 M , 
~~ Ay (Xp + Ху + Ха) + ñ 21 + X2 + X ) 





ird 
`x = 
Ay + > 
- 1 |(A61+42)X2+ (61+42)X4 + A1%3 + 42%; 
= X Ä ш Ls 
& 1 + À» 
AyXz + Ao, e 
- * | Mo + Ма + ы. сыч. (24) 
À | + À 2 
Similarly 
| | à Yz + À> Y 
ЕЕ e (30) 
Ay + А? 
= Е 2 Ат: +A2é M 
Pd [hene = n (31) 
Ay + Ao 


0.4 COORDINATE TRANSFORMATION 


Now that the local coordinate system is formed and properly 
located at the centroid of the element, the global coordinates of the 
corner points (X, Y, Z) are transformed to local coordinates (£, n, L) 


through the components of the reference vectors of the local coordinate 


sustem as follows: 


> ë 
nz. ШЕП nz. A Т X © 
Nz, n2, No, ET = 1 (52) 
Nn, Пу n; if ae ( 
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The corner points are projected into the plane of the quadrilateral 
element by setting the & components to zero. The original diagonal 
vectors, Tiz and T24, will be on opposite sides of the resulting plane. The 
plane quadrilateral element is now completely defined. The program will 
sweep through all of the input elements using the assigned location 


numbers, and repeat this process for each element. 


Figure 11. Forming the plane quadrilateral element 


о.о MOMENTS OF INERTIA 


The calculation of the moments of inertia for each element are 
performed for use in the computation of the velocity coefficients using 
the quadrupole method. Any calculus text will give the moments of 


inertia of a planar section with a constant unit density about the origin 
to be: 


= Th oe 2 u т z) 
I 7 8 £ 0 cd f| (55 


A 


Іш 7 H nó de dí, (34) 
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For the triangular region defined by the corner points 2, 3, and 4, 


es г ? ^ 42 
x = E A Es ea + 25761 (36) 
A / ? г 2 ғ ^ Z 
yy = Fp lUn2* mi * img * mí P *Cní * m 7] (37) 
А à ë x vu! 
Lay = 12 [( ë, + ë; )( Kerle ie) ШІ 1.50) 


ОЕ АИИ 


Similar equations canbe generated for the triangular region defined by 


the corner points 1, 2, and 4. The moment of inertia for the entire 
quadrilateral is the sum of the corresponding expressions for each of the 


triangles. Theresulting equations are: 





^ 
ға : ^ A V^ 
bx 7 42 E eal ee Es JE] ЕЕ! 
- ге ^ A x 2 
та УЕ еа EIER ye] 
Іш = 5 102 + Ne E ио 
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Ixy = Tg 10629 680 теж тз) +0659 ЁС та) (40) 
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6.0 THE MATRIX OF INFLUENCE COEFFICIENTS 





with the quadrilaterals completely formed, the next step is to 
calculate the velocities induced by the elements at the centroids of all 
the other elements. The total number of elements forming the surface 
will be represented bu N. Let the source element be the (j)th element, and 
the element for which the velocity components are to be calculated at the 
centroid is the (i)th element. It does not matter how the (i)th elements 
are arranged in relation to each other as the sequence progresses. 
However, the sequence must be consistent as the calculations proceed 
from one source element to another. This program sweeps through the 
(i)th elements in the order of their location numbers, m and n. For each 
consecutive n line, the elements are swept in order of increasing m 


numbers. 


Theresult of the induced velocity calculations for the elements 
with unit source densities is an N by N square matrix of the values of 
induced velocities at each element, known also as the "matrix of 
influence coefficients." The XYZ potential flow program calculates the 
coefficients column by column, while the Hess program calculates them 
row by row. The advantage of one over the other depends on the method 
of later solving the matrix for the source densities. In calculating the 
influence coefficients, twenty-five quantities which describe the 
geometry of the source element are required to adequately calculate the 
induced velocity at the centroid of the (i)th element. These quantities 
Include the coordinates of the centroid in the global coordinate system, 


the elements of the coordinate transformation matrix, the local 
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coordinates of the corner points, the maximum diagonal, the area, and the 
secondmoments of the quadrilateral element. Additionally, the Hess 
program uses the coordinates of the null point, making a total of 


twenty-eight quantities for that method (Hess and Smith 1962). 


When calculating row by row, the first (i)th element is selected, 
containing the "null" point, and the influence coefficients are computed 
for all of the (j)th elements in sequence before proceeding to the (i+!)th 
element. This procedurerequires the twenty-five quantities for each 
(th element to be available for calculation of the influence coefficient. 
Because each of the N (Dth elements is used N times with this procedure, 
calculating the geometric quantities or retrieving the values from low 
speed memory would be very time consuming, since the calculations or 
memory access would need to be performed N? times. Therefore, it is 
more practical to have the values available in high speed memory, 
although large matrices may exceed the storage capacity of high speed 
memoru, imposing a limit onthe number of elements which can be used. 
The advantage to the row-by-row calculation is that solution of the 
resulting matrix by the Gauss-Seidel reduction method does not require 
transposing the matrix, which would be another time consuming process 
(Hess and Smith 1962). 


Another alternative is calculation of the influence coefficients 
column by column. This method calculates the influence coefficients bu 
sweeping all the (i)th elements for each (j)th element before proceeding 
to the (j+1)th element. Therefore, the twenty-five geometric quantities 


are retrieved from low speed storage only once for each (j)th element, 
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for a total of N times. This procedure is not limited by the capacity of 
high speed memory, and calculation of the coefficient matrix is 
approximately 10% faster than the row-by-row method (Dawson and Dean 
1972). This is the calculation method used by the XYZ Potential Flow 


Program. 


An influence coefficient represents the combined effects on one 
element of the velocity potentials of all the other elements comprising 
the body surface. For the quadrilateral element with a unit source density 
in the xy-plane, from equation (13), the potential at point P (x, u, z) due 


to the element is 


1 dé In 
ф = f т JP ——— (42) 
Т (x - EY + Cy - т) +2 





ғ 





The integrand, 1/r, canbe expanded іп a series about the origin in 
powers of £ and ņ. Each term of the series will contain the product of 
some powers of è and ņ with a corresponding derivative of 1/rg, where rg 


is the distance of the field point P from the quadrilateral origin. 





and let 





Then the series expansion through the second order term is (Hess and 
Smith 1962): 
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P= AW~ (MyWy + MW) + 1/20 yu way * 2l Wy tlg ws... (43) 


The subscripts, x and y, used with w represent the respective partial 
derivatives. This series represents the multipole expansion of the 
velocity potential, since each term canbe interpreted as a point 
singularity of a particular order. The first term is the potential at point 


P due to a point source of strength A located at the origin. The second 


term is the sum of two dipoles of strengths M, and M,, located at the 


y 
origin, oriented along the x and y axis respectively. The choiceof the 


centroid of the quadrilateral as the origin of the local coordinate system 


causes the first moments, М, апа My to be zero. Therefore, the dipole 


terms dissappear, and are not dealt with anywhere in the program. The 


third term is the sum of three quadrupoles of strengths lyy, ЕТ апа ПТ 


located at the origin. Kellogg (1929) shows that this second order 
approximation is absolutely and uniformly convergent, and Hess and Smith 
(1962) show that convergenceis rapid enough with an increase in rg that 
certain further approximations may be made without significant error at 


large distances rg from the source quadrilateral. 


Hess and Smith (1962) presented a comparisonof velocities 
calculated using the exact formulas, a simple point source, and a second 
order approximation. The comparisons were based on the ratio of the 
distance rg, between the centroid of the source quadrilateral and the field 
point P, to the length of the maximum dimension t, of the source 
quadrilateral, typically the maximum diagonal. The non-dimensional ratio 


is then rg/t. Theresults show that sufficient accuracyis maintained 
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while using a simple point source at ratios of (rg/t) 2 4, using the second 
order source and quadrupole solution for the range 2.45 s (rg/t) « 4, and 
using the exact solution for ratios of (rg/t) « 2.45. Іп any case, the error 
goes to infinity as the field point approaches the edge of the quadrilateral 
where calculations indicate an infinite velocity. The XYZ Potential Flow 
Program uses a monopole source for (rg/t) > 4, the source - quadrupole 
formulae for 2 « (rg/t) < 4, and the exact formulae for (rg/t) x 2. Hess 
and Smith (1962) reported a maximum error of 0.001 in approximating anu 


velocity component using the above criteria. 
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7.0 DERIVATION OF THE EXACT SOURCE PANEL INTEGRATION 





From equations (2) and (42), the components of the velocity at the 


field point P(x, y, Z) due to the source quadrilateral are: 





wer 5 Eum (44) 
ox OY r3 
д 
d = ë 
У E - JP {y= de oT (45) 
OY r3 
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(46) 
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Equations (44), (45) and (46) are evaluated by expressing each of 
the integrals as the sum of four terms, each term representing the effect 
of one side of the quadrilateral (Hess and Smith 1962). This method can 
also be generalized for polygonal elements with any number of sides. The 
potential function for each side of the quadrilateral is the combined 
potentials of semi-infinite strips whose boundaries are the side of the 
quadrilateral and two semi-infinite lines parallel to either the x ory 
axis. When observed from the domain, and the sides are traversed in a 
clockwise direction, the source strip onthe right will have a source 
density of o = +1/2 and the sourcestrip onthe left will have a source 

_ density of o = -1/2 as shown in figure (12). When the sides are 
recombined to form the quadrilateral, the source densities outside the 
quadrilateral cancel each other, and the source densities within the 


quadrilateral combine to form a source density of o=+1. This will be 
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true for a planar element with any number of sides and in any relative 


orientation within the plane. 
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Figure 12. Fundamental potential functions for sides of a quadrilateral (Hess & Smith 1962) 


7.1 THE Y VELOCITY COMPONENT 


From equation (45), the velocity component Vy is found by summing 


the four terms representing the contributions of the sides of the 
quadrilateral. For the side from point (È, n) to point (£,, n2), the 
contribution is expressed as the integral over the area of the 
semi-infinite strips with the source densities of o = +1/2 and o = -1/2 


rather than the unit source density of equation (45). 
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Integrating with respect to т: 
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The terms evaluated at T| = +® and T = -co cancel, and the terms evaluated 


at T) = 112 add to obtain the following expression: 


t. 7 d£ 
V _ Bu pá 2 d£ | 
Bie 2 CTT КОО 2222 
2} Erg (us zi 
y) 8 ( š dé 4 - 
КЕЕ а = (4а) 
e| 


Equation (48) is changed to a function of arclength s by the relation 





e INS 
de EE Zu ee А. 
nm — í — (49) 
Р 5 2 ү di^ 





where dj? is the length of the side of the quadrilateral from (È m) to 


(82, 112) as shown in figure (13). 
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Figure 13. The potential due to a finite line source (Hess & Smith 1962) 


Substituting equation (49) into equation (48) 


dis 
>, > IE 
T $2 — А [ 
VI = 


sie dis J 


(50) 


ds 
r 


From figure (13), it canbe seen that, in terms of arclength s, the distance 


r from point P to any point onthe line from point I to point 2 is given by 





£2 — £i f К 
Баса | = (91) 
012 o 4114 а (5-5115 


42 





Evaluating the integral 


Vue = 











II 
=3 
C) 

ag 


The quantities гу, г», $, and s used in equation (52) are as shown in 
figure (13). Equation (52) is singular when r,-7 s, which occurs when the 
field point P is located anywhere along the line defined by the side of the 


quadrilateral. This singularity may be removed by using the law of 


cosines (Hess and Smith 1962). 





Figure 14. The law of cosines 
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From equation (52) 


e > 
| 22- 61 ro(i - cos fo) 
Y E O — o (cz 
Yie da 097.020 в.) (53) 


where B, and B; are the interior angles shown in figure (13). Applying the 
law of cosines to figure (13) 
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From equations (54) and (55): 
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Substituting equation (56) into equation (53) yields the final form of the 
exact equation of the y component of the velocity induced by the side of 


the quadrilateral frompoint 1 to point 2: 


= "+ Pe + dio 
i MUNIRI 
= di» fp + Pe dp 





Equation (57) is applied to the remaining sides of the quadrilateral 
simply by substituting the appropriate point numbers for the corner 
points of each side. The total contribution of the quadrilateral to the y 
component of the velocity is the sum of the four terms representing the 
contributions of each of the sides. The y component of the velocity at 


the field point P is now given bu: 











" pm | +2+ Û2 &2-& Го + Ге + Ооз 
у = O teas + BE lag e e 
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7.2 THE X VELOCITY COMPONENT 


A similar derivation process is used to produce the equation for the қ 
component of the velocity induced by the side of the quadrilateral from 
point 1 to point 2. The semi-infinite source strips are constructed 


parallel to the x axis, and the order of integration is reversed. 
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The x component of the velocity at the field point P due to the 


quadrilateral is given by: 


o. Fi €f» 02 3-1 r2 + Fz + daz 
Ww = —— n  lpq —— ———- in ا‎ 
012 Fi +P2- 02 025 Ze er 

(59) 
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7.3 THE Z VELOCITY COMPONENT 


The zcomponent of the velocity at the field point P due to the 
quadrilateral is obtained in a similar fashion, using semi-infinite source 
strips, this time parallel to the y axis. From equation (46), the 
fundamental velocity potential of the semi infinite source strips is 
integrated in a manner similar to that used to obtain equation (47), and 
the z component of the velocity due to the side from (EZ, 711) (0 (82, 12) і5 


given by 
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_ Performing the integration with respect to T, the integral 
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fits the integral form 


ағ - | F in- 3f dF 
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Then, from equation (60) 
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Again, the terms evaluated at +o and - cancel and the terms evaluated at 
N12 add to obtain the following expression: 


ae = » 
V A f E (U = Iis) d£ б (62) 
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Without a convenient substitution with which to integrate equation (62), 
the integration is performed directly. Recognizing that along the line 
defined by the side of the quadrilateral from (&, тр) їо ($2, 12), Тә тау 


be expressed as a function of 2; 


Tis = Miz Ë + Di (63) 


te 


where the slope of the side, my2 is given by 


"4-1 ` 
| te'h (64) 
MES уса a | 
e I 


and Dy may be determined knowing that niz = ny when È = è}. 


2 (65) 


So LU - Mié - Mo) 12 DR 
ao 9316 (00) 
т [0 = ё) +2 10 2) (Y Ma о) + 2° j Е 


Define the quantities 


Giz = Y = Da > M12* (67) 
и=х-ё (68) 

- Then du = -dé (69) 
U - Dio - My2É = 20 + фә (70) 





By a change of variable, equation (66) is now expressed as a function of u: 


dc (Mou + Qia) du 
м, „=2[ 


Ат Z 51722 5 (71) 
E). [иё + 2 Ци + (туи + фе) 22) 
DEG] 





X-E4 
i (Mis + Que) du 
20% SS a a le 79) 
ce [u 2 lim ¿+19 us 21M q u + qz + 2] ° 
лей 
which fits the form of 
| (Lu + M) du ae 
(Aue + 2Bu + Chal (аце + 260 + с) ros 
where 
L=mp Hauke 
A=] a = (mp? * 1) 
B = 0 р = M2912 
С = 22 Son 


From Hardy (1944), this integral form may be integrated by the 


substitution 
Ht +Q 





2 pra 1 
Ша (74 
where u and v satisfy 
аро + 6(р+о) + с= 0 (75) 
Аро + В(р+ о) +С= 0 (76) 
- and are the roots of the equation 
(аВ - DA) &¢ - (cA- aC) & + (bC - cB)= 0 (77) 
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Substituting the appropriate values into equation (75), the roots of the 


quadratic equation are 





I die (78) 
p m 
12 
Nh. Ze 
Р (79) 
12 


It canbe verified that these values satisfy equations (75) and (76). 


Substituting equations (78) and (79) into equation (74) 


3 
Mee” Ust 











Чә Па (50) 
Шіс == = me 
t+ | 
qt misz? 
IONS 
du - ale Me dt (en 
(t + 178 


Bu substitution and a change of variable, equation (72) becomes a function 
of the parameter t. After simplification, the integral now fits the form 


of 
en O | 

n^ " - - ; a 

N) (oct? ts 8) Oe, 


` where 
K = -(mi2°qi2z3 + 2mi23qi23z2 * q2 m12) 


«5 424% тр742722 
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В = 2424 + пг242222 
B= i? + qi 2° 


8 = 2024 + 12424 + 211124412222 + 22424 + 12242222 


Equation (82) canbe rationalized by the substitution 


Шы и 
. е uidi 


from which it canbe shown that 


„2 
° = et 
( 1 = ER (Ё 
$ ta 
dt - | (1 = JEK | gy (55) 


Substituting equations (84) and (85) into equation (82) and simplifying 


yields the integral in terms of the parameter v: 





(ытты: (ас (86) 
which fits the form 
] Ко i E tar! (87) 
where а2 = В 
b? z («$ - B?) 
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Performing the integration 





Ed ЕУ $5 = K E xó - р 
lp Ic tan ^ E (88) 


From equations (80), (83), and the expressions for œ, $, 8, and 3 from 
equation (82), and after a considerable amount of algebraic manipulation 


and simplification, equation (88) becomes 











Fis 22 - Qiu ч 
=- — An tan! | E i | (22) 
: riya fA + 28) - (u Ë Т, U M" 


2 jU*« iu - Tu4* « z? 
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Finally, applying these results to equation (71) 


X-T 
E f E ШЕТТЕРІН 
NS Е.Ы 
m х-®| гео (тооч до) +25 





= -tan | ===> 
| 2 4и2+ (0 - т) + 22 


tan" | Пиз - ёт) + — 
= (а م‎ ————_ 
| ZACK = & 12+ (0 - Ту) + ze 





¬ 


E 


ad Mia Ll - £2)* « 22) (y - mos - ёз) 
EH EE s Ts s 
ч ео)» a2: 





Recall that when x =¢,, y=, and when x = £,, y = >. Then, for the sake 
of amore compact equation, define the following quantities: 
2 = (х - 61) + 22 22 = (к - 62) + 22 


hy = (y - nx - Ey) ha = (y - n2X« - 22) 


The quantities r, and r, are as shown in figure (12), where 





ЕЕ ИЕ E 





Substituting these quantities into equation (91) yields the exact 2 
component of velocity due to the side from point (€,,7 1) to (£2,72) in the 


form used by the XYZ Potential Flow Program: 


E [ Mets = А 4 Mio Ez - Ho | 
B. - tan _ -tan a 


2 {3 
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The total z component of the velocity at the field point P(x, u, z) due to 


the quadrilateral element is the sum of the four sides: 


E 


V 


rJ 


E Мно Вт = hi p 
- Тап - tan 


2 
k 


E 


L 


+ tan” 


Таран ша ГІН | 


| 


+ tan” 


ШЕ 

[ж] {е 

ш | | ر‎ | : 
n pen | 
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8.0 APPROXIMATIONS OF THE INDUCED VELOCITY 
8.1 QUADRUPOLE METHOD 


As previously mentioned, as the ratio of (rg/t) exceeds the value of 
2, then certain approximations may be made which greatly reduce the 
calculation effort otherwise required by the exact method. In the range 
of 2 « (rg/t) < 4, the XYZ Potential Flow Program uses the second order 
approximation of the potential described by equation (43). With the 
origin at the centroid of the quadrilateral, the first moments are Zero, 


and the secondorder approximation is 


P= Aw+ (1/2)(1 „u Wu, + 21 Way 4 yy uy? (94) 


where the first term is a point source of strength A, the secondterm is 


composedof three quadrupoles of strengths lyy. lyy. and lyy located at 


y yy 
the local origin, and the subscripts onw indicate the partial derivatives 


of w with respect to those variables as before. The quantity A is the area 


of the element, and the terms | , and ІШ are the respective moments 


KK ГТ 
of inertia of the source element given by equations (39), (40), and (41). 
To obtain the velocity components at the field point, equation (94) is 


- differentiated with respect to the coordinate directions giving: 


99 





LO 
сП 
— 


1 
Dex ot Lx Wyxy $ > luy Wxyy | ( 


hi 


ip ES -[ Aw, + 


) 


a t 1 * 
= [ ANY y * 7 I Wyxyt Ixy ШТІ u Tyy Wyyn | (96 


nh] 


P Mm [su YYxuz + + [uu Wuyz | (97) 


кә|— 


Vz= шерте” = - Гду, + 


Recalling that w = —— d 


the derivatives of w, as expressed by Hess and Smith (1962) and as used in 


the XYZPF program, are 


W= -X g^ | 
Y 7 ZU m > (99) 
Wim 2 107 | 


= 


t; = =p 2 ж. 
Mr = © + 10% = Po d 
D y = OY ET, i 


A IS аңар 
“ыш 2 2541 


324 5 1044007! > (100) 


~” 

~ 
«C 
С 
ic 

Hl 


Cal 


Z pro? 


- Where 
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8.2 MONOPOLE METHOD 


When the ratio of (ro/t) is greater than 4, then the quadrilateral 
may be approximated by a simple source corresponding to the first term 
of equation (43). Then the velocity components at the field point due to 


the quadrilateral are given by 


op 
V. z — — = — AW, {111 
Ж Зу x ( | = | é 
a 7 
ل‎ = ——— = ВИ, 
/ y a у Aw 9 I 4 8 < J 
аф 
ME == =-AW, (105) 
dZ 


where the partial derivatives of w are those given in equation (99). 
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9.0 SOLVING THE MATRIX EQUATIONFOR SOURCE DENSITY 
9.1 JACOBI'S ITERATIVE METHOD 


From equation (20), the matrix equation may be solved for the 
constant source density O; for each element which satisfies the boundary 


condition equation (11). Equation (20) suggests the use of Jacobi’s 


iterative method of matrix solution in the form 


: N 
imt) I S I ~ ‚ ТК 
er = Vi + 3 Lii i ' 1 = Із Т, ... y N і | 04) 
ve: 
lz 


where N is the number of elements composing the bodu surface, and m is 
the number of iterations completed. Apartial sum of equation (104) is 
computed for each of the ith elements before proceeding to the next jth 
element. The iteration is complete when the summation of equation 
(104) includes all of the jth elements. Because the values of the source 
densities at all of the elements are recomputed before any of them are 
used in the iteration, this method 15 also called the simultaneous 
displacement method (Ralston 1965). This is contrasted with the 


Gauss-Seidel iterative method used in the Douglas program. In the 


Gauss-Seidel method, as each new Oo; is computed, it is used immediately 


in the iteration processfor calculation of Si) . This is also known 25 


the successive displacement method and is expressed as 
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(m+1) Bn (m* 1) à (m) 

= i | 
О; = Wi + = 1,0 +> Li; S; (105) 
1-1 


ішін! 


Е, М 


Though the Gauss-Seidel iterative method is faster, the Jacobi iteration 
method was selected for use in the XYZPF program in order to be able to 
perform the iterations column by column, since the coefficient matrix is 
also computed column by column, and the matrix does not have to be 


transposed for solution. 


When the (m*1)th iteration is complete, the values of the source 
densities are compared with those of the (m)th iteration and the 
differences summed for all of the elements. The total difference between 
successive iterations is then comparedto a convergencecritería input by 
the user. If the difference is less than the convergencecriteria, then the 
matrix solution is complete and the values of the source densities are 
stored for later use in computing velocities and pressure coefficients. If 
the convergencecriteria is not met, then the iteration processis 
repeated. After every five iterations, if the convergencecriteria is still 
not met, then an extrapolation is attempted in order to accelerate the 
convergence. The XYZ Potential Flow Program uses a Richardson 
extrapolation method, a numerical procedure which uses two approximate 
results to obtain a third approximation which is closer to the exact 


_ solution (Ralston 1965). 
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9.2 RICHARDSON EXTRAPOLATION 


The Richardson extrapolation assumes that the iterative process is 
convergent. For the iterative solutions 50, 5, and 55, Where Sy is the 
most recent approximation and S» the oldest, the solution is convergent 
if 


on 23 T 





х < 1 (106) 


91 - 02 


While a Richardson-type extrapolation cantake may forms, the XYZPF 
program uses a Procedure developed from the following approximations 
(Dawson and Dean 1972). If there is only one dominant eigenvalue and a 
sufficient number of iterations have been completed, the iterative 


solutions may be approximated by 


So 5; + EA" 
S¡ 5 Sp + E AM) (107) 
S) E Sp + E AM 


ES m-i 
ES +EX 


where Sr is the true solution 


A is the eigenvalue 
E is the eigenfunction 


m is the number of completed iterations 


Define the linear combination which, from equation (107), may be 
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approximated as 

Або + (1 - AJSS&SE EXIT (AX € 1 - A) (108) 
The value of A may be chosensuch that 

(АЛ +1 - А)= 0 (109) 
Then, from equations (108) and (109) 


A Sp * (1 - A)S * Sr (110) 


where the expression on the left convergesto the exact solution. 


From equations (106) and (109) 


ТІ 
Solving for A, 


= эу ЧЕ 
Шо t 4 ^N 


Ән di* X5 Ü 
Since the value of A generally changes from element to element, a 


weighted average of A is used in the extrapolation, where 


М 
2, (8200 - 8103) (sign of DC) 
=] 


PI 
n 
pd 
mo 8 
f, 
— 


М 
> Di) 
1=1 


Equation (113) is recomputed after every fifth iteration. 1f the difference 
between the new value and the old value is less than 0.02, then the 
- solution is extrapolated. From equation (110), the extrapolated solution 
is 
8* - А90 + (1- А) 5; (114) 
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When there are two dominant eigenvalues, then the iterative solutions 


may be approximated by 


5. x 5, + E, Amm! т Е, A7! (115) 


where Sr is the true solution 


Ay and A) are the eigenvalues 
E, and E? are the eigenfunctions 


m is the number of completed iterations 


Define the linear combination which, from equation (115), may be 


approximated as 


B, So + B, T + (1 - В; - В) 52 
^ 9f + E SS [Вә Ay? Í В; Ay + (! = В; n В.) (16) 


+E, AM? [Bo A? + B, À> + (1 - B; - B,)] 


The values of B, and B? mau be determined for which the eigenvalues 2, 


and A» are roots of the quadratic equation 
В, 2+ В, А + (1 - В, - В,) = 0 (17) 
_ Then, from equation (116) 
Bz So + В, 51 + (1 - В, - В) 52% $; (118) 
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where the left side of the equation (118) convergesto the exact solution. 


Using equation (115) and eliminating terms containing E»: 


(50 = Sy) m А? (5) ксі So) - Е, ША AJA d 1) 
(84 - $2) - A2 ($2 - $3) ES AP 30.473230, - 1) (119) 
(S; 2: p» À» (S; = San - Е, Q URSI. AMA y = 1) 


Solving for A, 


w=. =, % Е 
DT I-A 


та О ES {сус 
ш M 2595 2) m 


Ti ) 
53-52) - %2(5.- 85) 152-54)- 52152-54) 


From equations (117) апа (120) 


ЫЗ ІП 2211151- 22)-(ӛз- 54)) 


(1211 
D x 
(94- 52)(54-254%52)-(54- 53)1184- 52)-ібз- 54))| 22. 
а a (122) 
D 
where D- (54 = 253 = 556 = 252 М 54) = (5, Е 52 т 23 + 54)? 


The weighted averages of B, and B» are used for the extrapolation as done 
with A in equation (113). If the sum of the absolute values of the 
weighted averages of B, and В, changes by less than 2%, then the 
extrapolation is performed. Then from equation (118), the extrapolated 


solution is 


— 
— 
n5 
C 
MÀ 


moe peat БІЗІ (| -Г, -В2 162 
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10.0 CALCULATION OF VELOCITIES AND PRESSURE COEFFICIENTS 


with the influence coefficients and the source densities 
determined, the calculation of velocities is arelatively simple matter. 
From equation (9), the total velocity is the sum of the freestream 
velocity and the disturbance velocity due to the body. The product of the 
source densities and the influence coefficients are summed for all of the 
elements, and then added to the freestream velocity to determine the 
total velocity at any point in the domain. Velocities onthe surface of the 
body are calculated at the null points only, as the boundary conditions are 
enforced only at the null point of each element, and velocities at other 
points in the element would produce significant error due to the method 
of approximation. The components of the velocity at the centroid of the 


ith element are 


Vi x Vooy I қ Cijx Si 
ILS 
N a 
Ey hi ETT را‎ 
Т = E оо y + 2. 8 1] y ej 
= 
М 
. - Мт И + er i5 
E Vous Lii. ] 


From equation (15), the velocity induced bu an element at its own null 
- point has a magnitude of 27r directed along the outward normal vector of 


the element. 
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Ata point off the surface of the body, the components of the 
velocity are determined just as if the point of interest was a null point of 
a single element. The total velocity at the field point is the sum of the 
freestream velocity and the contributions of each of the elements of the 
body surface. The contribution of each of the elements is determined by 
calculating the influence coefficient based onthe element geometry, and 
multiplying the result by the source density for the element. The total 


velocity at the field point may be expressed as 
N 
Vo = Yo E Са © (125) 


where p and q represent the field point and the source element 


respectively and the influence coefficient, 


| rn Ес 3 dr Drs 


AS discussed in section 6.0, the influence coefficient may be 





calculated by the exact method, or it may be approximated by the 
quadrupole or monopolemethod depending onthe ratio of the distance, ro, 
between the field point and the centroid of the source element to the 


maximum dimension, t, of the source element. 
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The magnitude of the velocity at either the on-body or off-body 


points is given by 





The pressure coefficient is calculated by using the result of 


equation (127) in equation (19), renumbered here for clarity. 








66 





i 


11.0 STREAMLINE CALCULATIONS 


For the calculation of streamlines off the surface of the Dody, a 
timestep procedure is performed by calculating the velocity at the 
starting point of the streamline from equation (125), and advancing the 
streamline one time increment by a fourth order Runge-Kutta integration 
to anew point (Ralston 1965). The timestep procedure is repeated, thus 


creating a streamline composedof finite segments. 


For the calculation of streamlines onthe surface of the body, the 
Streamline is started at aspecified point and quadrilateral number. The 
local velocity is calculated from equation (124), and the values of a 
stream function are computed for each corner point. The stream function 
is chosenso that it has a value of zeroat the last point onthe streamline 
In the quadrilateral. The side of the quadrilateral through which the 
streamline exits is determined, and coordinates of the point onthe side 
which has a stream function value of zero are computed. The direction of 
the streamline is verified by comparing it with the known direction of 
positive velocity. The next quadrilateral through which the streamline 
passes is determined by calculating the proximity of the new 
quadrilateral to the most recent point onthe streamline. Acircular area 
is computed which encloses the new quadrilateral with an additional 10% 
margin. If the last point of the streamline falls outside the circle, then 
the quadrilateral is discarded and a new one selected until the streamline 


IS adjacent to the new quadrilateral. 
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Figure 13. Calculation of on-body streamlines 
This procedure is repeated along the surface of the body until all of 


the surface elements have been tested. Theresult is a streamline 


composedof segments from one side of an element to another. 
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12.0 DEVELOPMENT OF HIGHER ORDER PANEL METHODS 


The XYZ Potential Flow program assumes a constant element 
source panel as described in Section 3.3. Extensive use of the constant 
element source panel method has shown that the primary diasadvantage of 
the method is that, in order to obtain a highly accurate solution, a large 
number of surface elements must be used to discretize the body surface. 
The method has been applied to problems of increasingly complex 
configurations (Hess 1977). By doing so, the size of the coefficient 
matrix is increased resulting in increased computer time and cost. 
Additional cost is accrued due to the manhours required to prepare the 
input. Therefore, while the constant element methods have provento be 
very successful, the cost has motivated the development of higher order 


methods. 


The higher order surface sinqularity methods discretize the body 
surface with curved elements having a variable source density, as 
comparedto the flat elements of constant source strength used in the 
basic method. Hess (1973) showed that the effect of a curved surface and 
the effect of a variable source density are of the same order of 
magnitude. Therefore, the two effects must be used together to provide a 
“consistent” solution. The consistent higher order panel method provides 
the increased accuracyand speed desired for three dimensional Neumann 


- problems (Hess 1979). 


Accordingto Hess (1979), the evolution of the higher order panel 
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method from the constant element method involved the derivation of new 
influence coefficients based on the integration of a variable source 
density over a curved element. Other portions of the method were 
unchanged. However, the development of the higher order velocity 


equations also required different programming logic. 


In examining the potential for development of the higher order 
methods, Hess (1979) noted that "a consistent approach always uses a 
source polynomial one degree less than the panel polynomial." Through an 
independent effort, Brebbia (1984) presented a higher order approach 
using the direct method to solve for a surface potential polynomial 
stating that the potential function must be of a degree at least equal to 
the degree of the polynomial describing the element. Knowing that the 
velocity function is the derivative of the potential function, these two 
observations agree. Asaresult of his derivations, Hess showed that the 
solution of a flat element with a constant sourcerequires one integral, a 
paraboloidal panel with a linearly varying source density requires six 
integrals, and a cubic element with a quadratic source density requires 
twenty-three integrals. Development of higher order methods nas 
focused on the paraboloidal element with the linearly varying surface, as 
solutions of higher order than that offered little benefit for the amount 
of effort required to produce a working program (Hess 1979). Hess 
(1979) and Eriksson (1983) have independently developed programs for 
_ three dimensional higher order panel methods. The higher order Hess 
program evolvedfrom the constant element program which he developed 


in the early 1960s, while Eriksson developed a new program based on the 
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work of Johnsonand Rubbert (1975). Continued work in the near future is 
expected to deal primarily with refinement of the paraboloidal element 


with a linearly varying source (Eriksson 1983). 


In order to alleviate the burden of preparing the input, a geometry 
package for input data generation has been developed which is 
incorporated into the Hess higher order panel program. This allows the 
user to enter relatively few points to describe the body. The geometry 
package enhances the surface representation by distributing additional 
points onthe surface based on one of many algorithms or recurring 


geometries (Halsey 1978). 


As the state of the art in fluid dynamics has progressed, the XYZ 
Potential Flow program has seen increasingly complex applications 
requiring a great deal of effort in preparing the input, and requiring long 
computer run times. Hess (1979) reported the use of the Hess constant 
element program for a configuration utilizing 7000 effective elements. 
Realizing that the computation time increases as the square of the number 
of elements, it is easy to see the motivation for developing the higher 
order panel methods. Though modern computers offer storage capacities 
which can handle most applications of the constant element method, the 
higher order panel methods can provide equal accuracyfor much less user 
effort. While the constant element method is still a versatile tool, 

- future generations of the surface singularity methods will be able to 


handle the more complex applications being demanded in fluid dynamics. 


71 





| 


13.0 VELOCITY CALCULATIONS FOR A TRIAXIAL ELLIPSOID 


As the only true body for which an analytical solution exists, a 
triaxial ellipsoid was selected for the sample calculations in order to 
compare calculated results with the analytical solution. Hess has made 
use of the triaxial ellipsoid throughout his works in developing both the 
constant element method and the higher order panel method. Therefore, 
the XYZPF program will be compared with existing results of the Hess 


method (Hess 1979). 


The triaxial ellipsoid utilized for the calculations has semiaxes 
dimensions of 1, 2, and 0.5 in the x, y, and z directions respectively. The 
surface was discretized by selecting fixed intervals of 0.1 in the y 
direction, and fourteen equal divisions of the 90° sector in the x-z plane. 
The values of x and z were then solved in terms cf y and an angle 8. This 
method yielded 280 elements in the first octant for a total of 2240 
effective elements after employing symmetry. AFORTRAN program was 
used to generate the corner points and to prepare the input file for later 


use by the XYZPF program. 


Figures (16) and (17) show excellent correlation with the analytical 
solution and little difference from the Hess solution using 4320 effective 
elements. Theuse of the centroid as the control point is an 
approximation used to simplify the multipole expansion of the potential 
about the origin of the local coordinate system. This approximation is 


valid for most elements. However, for elements which are long and thin, 
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the physical difference between the centroid location and the null point 
location is significant, and use of the centroid can produce significant 
error as may be observed in figure (16) when the value of y approaches 


2.0. 


Recent calculations on the same body (Hess 1979) showed that 
results of at least equal accuracy could be obtained using only 480 
effective elements using the higher order panel method. These results are 
a significant demonstration of the value of the higher order panel method. 
Using the higher order panel method rather than the constant element 
method, the user nas the option of obtaining equal accuracy with cruder 
discretization or higher accuracyfor the same discretization effort. 
While the results of the triaxial ellipsoid show relatively little 
improvement in accuracy, the most significant advantages of the higher 
order panel method are evident for a body with concave regions (Hess 


1977). 
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ANALYTIC SOLUTION 
o HESS PROGRAM (4320 ELEMENTS) 
С]  xYZPF PROGRAM (2240 ELEMENTS) 





Figure 16. Comparison of analytic and calculated velocity distributions on an ellipsoid with 


axes ratios 1:2:0.5. Velocities in the xz-plane. (from Hess and Smith 1962) 
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ANALYTIC SOLUTION 





о HESS PROGRAM (4320 ELEMENTS) 
EJ  XYZPF PROGRAM (2240 ELEMENTS) 





Figure 17. Comparison of analytic and calculated velocity distributions on an ellipsoid with 


axes ratios 1:2:0.5. Velocities in the yz-plane. (frora Hess and Sroith 1962) 
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14.0 СО | М MA 


The objectives of this paper were (1) to describe the details of the 
approximation of an arbitrary three-dimensional body using quadrilateral 
elements, and (2) to provide a detailed derivation of the exact source 


panel integrations. Both of these objectives were met. 


The method of surface discretization and source panel geometry is 
easily described using basic principles of geometry and vector algebra. 
By using quadrilateral surface elements, many surfaces can be discretized 
ina very straight forward logical fashion. The user can frequently 
visualize the contour lines of the surface which may be used to form the 
quadrilaterals, with some help from an intuitive approachto the fluid 
dynamics problem. The method of forming the planar quadrilateral 
element in the XYZPF Program differs slightly from the method presented 
by Hess and Smith (1962). The differences lie in the formation of the 
local coordinate system and the use of the centroid rather than the null 
point as the control point for applying the boundary conditions. The 
method of forming the local coordinate system has no effect onthe 
potential flow calculations as long as one of the coordinate vectors is 
the outer normal to the planar element. The use of the centroid as the 
control point is an approximation used to simplify the multipole 
expansion of the potential about the origin of the local coordinate 
system. This approximation is valid for most elements. However, for 
elements which are long and thin, the physical difference between the 


centroid location and the null point location is significant, and use of the 
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centroid can produce significant error as may be observed in figure (16) 


when the value of y approaches 2.0. 


A detailed derivation of the exact source panel integration has not 
previously appeared in literature, though the results are summarized by 
Hess and Smith (1962). The derivations presented in this paper verify the 
equations presented by Hess and Smith (1962), and the equations used in 
the XYZPF Program. The integral expressions for the velocity 
components were evaluated exactly with no assumptions or 
approximations used in the course of the integrations. Since the method 
of integration reduces the surface integral to a line integral around each 
of the sides of the element, the integration method canbe generalized for 
a planar element with any number of sides, though the surface 
discretization used in the XYZPF Program uses only quadrilateral 


elements. 


The calculation of potential flow about arbitrary three dimensional 
bodies is an engineering tool which is basic to design involving fluid 
dynamics. The XYZPF Program is a useful tool which has proven its value 
over the past 14 years. However, the increasing demands placed on this 
method are exposing the errors of the approximation as evident in the 
sample calculations presented in this paper. The requirement for 
increased accuracyhas motivated the development of the higher order 


panel methods. Some of the limitations imposed onthe XYZPF Program 


were due to computer memory and speed limitations. Advancesin 


computer performance may allow future investigators to eliminate some 
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of the simplifying approximations used in the XYZPF Program, allowing 
increased accuracy without violating computer limitations. Some 
modifications might include the use of the null point as the control point 
rather than the panel centroid (as is used in the Hess program), or 
extending the range in which the exact velocity calculations are 
performed. The gains in accuracyby modifying the "constant element 
method” are limited by the basic approximations of the planar element 
and the constant source density for eachelement. Significant gains are 
most evident in the higher-order panel methods. This author concurs with 
Eriksson (1983) in expecting advances in the surface singularity methods 
to focus onthe “development and refinement” of the higher order pane! 


methods. 
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Appendix | - XYZPF Section PF! 
Reads Input and Computes Era 
Parameters 


Appendix П - XYZPF Section PF2 . 
Computes Influence Matrix Coefficients 


Appendix Ill - XYZPF SectionPF3. 
Solves Matrix Equation for Source Densitu 


Appendix IV - XYZPF Section PF4. 
Computes Velocities and Pressure Coeff icients 
for On-Body Points 


Appendix V - XYZPF Section PF5 . 
Computes Velocities and Pressure Coeff pie 
for Off-Body Points 

Appendix V] - XYZPF Section PF6. 
Computes Velocities and Pressure Coefficients 
for Of f-Body Streamlines 

Appendix VII - XYZPF Section PF 7. 
Computes Velocities and Pressure Coefficients 
for On-Body Streamlines 


_Appendix Vill - Triaxial Ellipsoid Input File . 


Appendix IX - Triaxial Ellipsoid Output File . 
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APPENDIX І - XWZPF SECTION PF1 





C) 652 o O 


Ье 


10 
20 
50 


ж”, 
sn 


E 


PROGRAM FFFIXINFUT= 128, OUTFUT= 128, TAFES= INPUT, TAPEG=ÖUTFUT , TRFEOS, 
TRFES=TAPEQS , TAPEQ4 , TRPES=TRPEOS , TAPESU= 128 5 


X¥Y2 POTENTIAL FLOW PROGRAN VERSION 4 SECTION 1 
READS INPUT AND COMPUTES QURORILATERAL PARAMETERS 


FOR INFORMATION CONTACT 
BILL CHENG OR JANET DEAN 
NUMERICAL FLUIO OYNAMICS BRANCH CODE 1643 
HRURL SHIP RESEARCH ANO DEVELOPMENT CENTER 
EETHESDA, MARYLAND 20054 


DIMENSION INDEXCS, 3>,G¢9,6), F699, C209), FQ) XPCQ), YPC) ZP CG) 
,HSKC1005, NS¢240>, PROBC 15>, DHCSSO 

COMMON XCS00 >, (GOD 9, 26800 25,1D£4 1,215, BC2500, TC45000, KP C 1002 

EDUIURLENCE — CCZCTO,FC 1» ) 

EQUI URLEHCE£WEX 12, KF€ 132, CHSCTO 12, HEK C120, CHE C20 DD, NP. >, CRC ZO 3, 

NSP), CHEC2D032, NEP, CHEC204 2, NSE », CHE C2052, MLXO, CHE C206 5, HL 


و“ 


ИЕ ЖООГО Кра, |с, (с?й? [®Н12,(И40211›,К›,сызєт13), 
4 EPS) CWECZOR), IPS У, СЕСТРУ, 1РЕ 3, HSC2 120, XD, OIS C2 ED, VE 2, 


45:2 192,212 


EQUIIURLENCE. ($12,232, (84, 415 
INTEGER P,P1,P2,P3,F4,PC,PS5,PB,P?,PB FO 
WRITE (6,5) 
FURNATC49H 1x42 POTENTIAL FLOW PROGRAM SECTION 1, VERSION 4 3 
FORMAT (11, 1549) 
FORMAT (1X, 1514) 
ШІБІНІ СТХ/217,БЕ12,52 
FORMAT (1Х,5Е12.9) 
А. FERD IH CONTROL PRRRHETERS 
HSC 2s =<. 
K1zn 
Il 1=0 
02=0 
05=0 
ІГІ4- 
|O5=0 
|DE=0 
[ÛU ?=Û 
МАХМ= 20 
NAXN=40 
МАХНОЕ=6 50 
MAXPC=S00 
ICNTRL= 1 
EGF SO=0, 
READ «5 , 102J, «PROBCI, I2 1, 152 
E EUrCo» .EQ. O.) GOTO 9 
UBITECE, 8) 
FERMATCSOHING TITLE CARO FOUND - PROGRAM AEORTED 2 
STOP 
CONTINUE 
J=Ü 
WRITE «8,10? J, <РРОЕСІ2,1-1,15» 
<Н-.0 
5Б-0. 
FORMAT CIPHARD. OF QUADS. =,14 ¿VOR WO. DF SECTIONS=, 


214/31H MAX. NO. OF ITERATIONS X FLON ,12,9H Y FLOW, 
313,9H 2 FLOW ‚19 


READ <5 ,119NQE, NSE, MIX, MIY, MIZ, ISM, EFS, IUCT, IPS, IPF, ISP 





1 ЕЕЕ Е ЕВИЯ ТАРЕ ХСЕНТЕЕ , ҮСЕНТЕН, ZLENTER 
11 ЕОВМАТ (614 ЕЁ. 5,814, 1X, 3F5.3) 
IF <EOF£S> EQ. 0.2 GO TO 19 
WRITECG, 18> 
18 FORNATC43HONCG PARAMETER CARO FOUND - PROGRAM REORTED ) 
STOF 
19 CONTINUE 
IF CIEDIT1.E0.1y LECNTRAL=0 
WRITE (6,21) NOE, NSE, MIX, MIY, MIZ 
БІРЕПЕТІНТ ( CONVERGERCE CRITERIA” FS. DO 
41 FORMAT C4HG M,7X,2HX1, 12Х, 2НХ2, 12%, ZHXS, 12X, HIS EM: 
ШОНУН, 12Х, en N,?X , 2HV1, 12X, 2Hv2, 12X, 2Hv2, 12X, 2HV4, 
PX НУР, 12X,2HVN, 12X, 2HFL, 12X, 3HCZ5/4H — P,?X, 2H2 1, 12X, 2HZ2, 12X, 
3 2H22. 12X, 2H24 , 12X, 2HZF , 12X, 2HZH, 12X , 4HCZ 1 | 10X, 3HC25/ 
ЖЕБІР СЕЛІ НУ 201, 125,252 12Х,2НХ3, 12Х,2НХ4, 12Х,2НАР, 12Х, 
| ZHXH, 12X, IHR, 12 3X, SHCZ4 /4H N, Px, 2H1, VER, 2ҢУ2, 12X, 2HY2, 12X, PHYS. 
IAL, НУР, 12x. ZHYN, 12X,ZHFL, 12%, SHCZS 24H m ӨНУІ 127-242) ter, 
2 2H23, 12X, 2Н24, 12X, 2HzP, oo. 2HZH, 12X, 4HC2Z 1 ‚10x, 246267) 
24 FORMAT CH , 11, 1H PLANES ПЕ SYMMETRY 3 
240 WRITE (5,24) [SH 
КОНЕЦ ТЕ (6,315 EPS 
zen IF CIFS.LE.U0» GO TE 200 
Pee HRITE (6,562 IFS, IFF 
26 FORMAT C4SHOHNEN SOURCE nue D jn ES ru TEU FOR А05 514 SH = 
1147» 
290 К=0 
ЫН!ТЕС6,ЗӘМЕР 
sa FORMAT CcaHü ISP = 150) 
ШЕ 572 IEBIT1 IEDITS IEDBIT4, ITHPE 
FRAMES СӨППІЕПІТІ -,ІЗ/АН ІЕШІТЗ - ІЗЖН IEDIT4 =,13/9H ITAFE =, 
1 13) 
ШҮЙЛЕКЕ 2C) ХГЕНТЕН,ҮСЕНТЕН, 2СЕНТЕНВ 
> FORMAT <¢ 10НОХСЕНТЕВ ecc MUN NEENIER T ES 7/UH ZECEHTER z.F5.2) 
ММ=П 
MN=Cı 
P=] 
Q=1.0 
DO 29i |-1,41 
ПП 291 d=1,71 
291 ІШС1,40-0 
4-0 
EAE RENOU FIRST PT. 
IERF=Ü 
[F CITRFE.EQ.1> GO TO 292 
EE BERL (9,405 XI, VIL, 21, NI, HL, HS, HE, UN 
IF (EOF (5). NE. 2. “OF. NE. UE. 0» G0 TO 2050 
WAITECG 45) ICHTEL М5 
45 FORMAT C11, ан SECTION , PE 
LIHEZÜ 
iu ТП 293 
2050 IECIERA.E0.09 бб ТО 2200 
WRITECE, 21005 
2100 FORMAT CSOHONG POINT CRAIC FOUND - PROGRAM ABORTED > 
STOP 
2200 IEFF =1 
ITRPE=1 
WRITECE, 2300) 
300 FORMATC4?HDERROR IN INPUT - POINT CAROS NOT ON INPUT FILE, 10X, 
1 SOHPROGRAM WILL CHANGE {TAFE TO 1 AND TRY TO REAT TAFESO > 
IF < ЕСЕ(52.Е0.0 2 URITECE, 2400) ү 2! 
2400 FORMAT 11HOEXTRA FLOW, ЗЕ [2.5 Эу, 2OHHILL NET EE CONFUTED > 


2 


сә 
on 





292 
40 


2000 


293 


sna 
301 


304 
305 
306 


РЕНП(50,402 ХІ,УІ,21,М1,М1,М6,МЕ,Ч 

FORMAT (3F12.9,414,F12.9> 

FOF SC=E0F (S03 

IF CEOFSO.NE.O. .OR. NS.LE.03 GO To 2450 

WAITESG,459 ICNTRL,NS 

LINE=0 

GO TO 293 

IF СЕНЕР ЕП 02 60 TO 2500 

WRITECG, 2100) 

STOP 

[ERR= 1 

ІТНРЕ-0 

WRITECE, 26009 

FORMATC"ERROR IN INPUT - POINT CARDS НОТ ОН TAFESO", 10%, 
("PROGRAM HILL CHANGE ITAFE TG O ANDO TRY TO READ INFUT FILE") 
GO 10 2000 

UNR=UN 

HSS=HS 

FC= 1 

IF «NE .EQ.O> бб Т0 2700 
НМІН-НІ 

HHRxX=HI 

NMIN=HI 

NMAX=H 1 

GU Т0 500 

МНІН-НІ 

MHAX=H ! 

ЧИТЕН 

НМАХ=М | 

GO TG 300 

IF CITAFE.EO. 13 GO TO 2237? 
БЕН 55 407 ХІ/ТІ,ГІ,НІ,ПІ,М5,ПМЕ,ЦМ 
ІЕ «ЕОЕ(52.Ей.0.2 б0 ТО 204 
5-0 

ШЕГІ, 

Is. 

21-0. 

BO TO 299 

БЕНГ 50,40) ХІ,УІ,ГІ,НІ,НІ,М5,МЕ,ШІН 
EOF SO=EOF CAO) 

IF CEOFSO  .NE.0> M5=0 
BESDPU 1 

ЖЕНЕ ТЕ Hes) GO Te 320 

IF СМЕ .ЕП. 0) GO TQ 304 

IH=NI 

МЕНШ 

Мізсін 

Ері SIGRE РТ. ІН РТ. RARAY 

ШЕ С ПХРЕЗІ-РС2 205,305,310 
МЕТТЕСБ, 305) МӘ,МІ,МІ 
FORMATCGQH ERROR IN INFUT = THERE ARE TOC MANY DATA POINTS IN SEC 


ITION ,19, 30H - POINTS BEGINNING HITH M =,14,5H N= ‚14, 


2 IH HILL EE IGHOREO > 


LINE=LINE+1 

р4=1 0144 1 

GO TO 295 

XACPCOZXI 

YPC = | 

LEPC I=? | 

fee ih be Hay AND. МІСЕ. Махмэ GO Td 315 
WAITECG, S112 MILNI 


Гу) 
сл 





LINE=LINE+1 
311 FOPMATC3SH ERROR IN INPUT = INVALID M,N INDICES —,10X, 
1 14HPGINT WITH M= ,14,5H N= ,14,17H WILL BE IGNORED > 
1DS=/D5+ 1 
EEE 
bo TO 295 
S15 (OME, NI SPC 
MMAX=NAXQCMMAX , MLD 
MM IN=NENOCMMIN, MI > 
NMAX=MAXOCNMAX, NI > 
HMIN=MINGCHMIN,NID 
GO TO 295 
330 IF «IEDIT1.EG. 1» Gü TG 294 
IF CLINE.LT.40)> GO Tu sss 
HRITEC&, 42) 
LINE=0 
50 70 294 
333 HRITECH,41) 
294 CONTINUE 
E. 00 LOOFS To SWEEF PT. ARRAY 
NI=NMIN 
ММ2=МИАХ-ММ М 
NNZ=NMAX-NHMIN 


IF € MODCHN2,22.E0.O .AND. MODCNN2,22.EQ.0)? GO TO 222 
HRITECB, S21? NES, MN, HE: NET, INHRI: 
LINE=LINE+ 1 
331 FORMRTCIBHOERROR - SECTIOM ,15,45H DOES NOT HAVE QUADS ARRANGED IN 


1 ELOCKS OF 4 ,0H — HMINs ,12,6H HHfRXs ,12,6H НМЇН= ‚12, 
26H NMAX= ‚122 
106=106+ 1 
322 ММ2=ММ2/2 
NH2zHN2 /2 


Ut 404 Nid=1, NHZ 


МІ-ММІН 
DO 402 Mi=4, MHZ 
NO= 1 


F. HAVE & CORNER PTS. BEEN GIVEN 
IT2IDCH 1, N TE TDCH 1+ 11 LOCH +2, NT EID М ЖОМ Я, МЕ 6 
СИЕНИТ УУСУН 
IFC IT.EQ.0 > GO TO 402 
IERE-O 
M2-M 1e 1 
М2=М1+1 =М1, М2 
Ой 400 МеМ1,М2 
GÜ TO (334,535,326,2337) МО 
324 PlsIDiH GN 2 
P2sIDéMe 1,M 9 
PZ=1EICH+ 1, N+1) 
Р4=104М ‚N+1) 
PS=IDCH+2 N > 
F&sIDZHZ, He 12 


P?sIDCHs 1, N42) 
FO=ICKM N42) 
PG-F 1 


IFXXXXP 12. HE. XCP22. OR VXP DX. NE VXPZ S. OR. ZXP 12, NE. ZXP22).— AND. 
1 (OXCP D. HE. XCP42. OR CP 12, NE VCP4 D. OR ZCP 12. HE. ZCP4 20? GO. TO. 240 
Fü- IDCH*Z N42) 
GO TO 340 
335 PisIDXM — , H4 12 
P2-1DXM N 9 
F3=1DCM+1,N > 





1 


P4=1DCN+1,N+1) 
PS=IDCH N-19) 
F&= IDCM+1,N-1) 
P7=IO(M+2,N > 
РЕ=10<М+2,М+12 
РӘ-Р1 
IFCCXCP 42. NE. X€P22. 0R VCP 12, NE V(P22.0R.24(F 129. NE.Z4P22).— AND. 
(XCP 12, NE. XCP42.0R VCP 12. NE VCF4D. 0R. ZCP 12, NE. ZCP4 222. GO TO 340 
PGzIDCHe Z2, N-1) 
Go TO 340 
Pi=1DcH+1,N > 
P2=IDCH+1,N+1) 
F3=IDCH М1) 
P4=10¢H N 9 
PS=IDKM+1,N+2) 
P&sIDCM ,‚N+2) 
P?-IDZH- 1, N15 
P&=IDcH-1,N 2 
PO=P 1 
IFCCXCP 13. NE. XCP22. OR SCP T2. NE VCP2). DR. Z£P 1). HE, Z£F220.— .RND. 
CERO ME XCPA 2. OR VCP 0. HE VCFA 2. OR. Z0P 12. HE. 2CP4 222 GO TO 340 
PüsIDCH- 1, N22 
GO TO 340 
Р1=10<М+1,+17 
Р2=104М №12 
Рэ=10<М N > 
Р4=10<М+1,М ) 
PS-IDCH- 1, N15 
F&sIDCH- 1, HN > 
P?SIDXM ,N-1) 
P&- IDCM 1, N-12 
Pg=P 1 
CORO ME XCPZ 2. OR VCP 12, NE ODO OR. 2010 NE. ZCF2D) AND. 
(OX£P 13. HE .X4F4 9. DR V(P 12. NE, VCP4 0. DR. ZCP T2, HE. Z(P4 2)? DD. TO. 34D 
Foz |DCH- 1, N- 12 


340 ІРСІЭ-Р1 


[FE =F 
ІРСЗЗ-РЗ 
ІР" Ӛз-Р4 
ІР(О>-РЕ 
(PCE }=Fë 
IPC? =P? 
IFCS FS 
[FeO 3=PQ 
GZ COMPUTE NORMAL VECTOR (XN YN, ZMO 
AS CPT) 
X2zXCtP4o-XCPZ2)5 
V IEVCPSO-YCP 12 
EN CEA CPZ 2 
¿ds 2 PD Z AP 1) 
22-21Р4о0-2ҚРоо 
XNzVZ2]1-V 1522 
ҮЧ=Х 1#72-Х2#*7 I 
2Н=х2*' 1-х1%\'2 
R-EDZCXH, VN, 2H? 
IF «R.GT. .OQo000000001» GO TO 345 
ЫБ!ТЕ(6,3432 
FORMATCS3H ERROR IN INPUT = ZERD AREA QUAD > 
LINE=LINE+1 
¡Q01=/01+1 
Rü=l, 





n 


342 


945 


J F-J FJ СӨС 


GO TO 251 
CONT [HE 
XNZXH/R 
YN=YN ¿E 
¿n=2N ¿Rh 
Al=.S#F 
COMFUTE CENTROIOD 
Mela Р5)-дР2) 
ИЕ Ce SP 2 OD 
AAA 
X= Y 122 -P22 1] 
lxx o | 
Bo SY 
В ї=5002 059, УЕ, 29) 
HZ-R-H1 
IT=1 
a TPO UPS + CHIEN CPS AZ FAP IIA), 
MOH CV CPZ FN CPS DHA ИКЕ ОН? жәнне 
ESC ZIP LCP OA CATRE PSI RZF ZIP I DI FH 3 
CaMFUTE SECOND AND THIRD VECTORS 
XIV NFZI-VIFZN 
ісе 1-2 FAH 
ыле I-X THUN 
б-1./202(Х4,Ү4,74) 
Х&=х4*®Ң 
MdczWy4*H 
24=24*R 
X= 29 ZN 
yV3=xH+Z4-X4*2H 
ESSE YEN 
REUMEUTE POINTS 1H QUAD SYSTEM 
00 947 1=1,9 
= 1Р1 
КЕ Елан JAC IHY IEY CL -YC Zat L-EN 
EIEC SIC IY EEVEE Jt Z4 FCECL I-T) 
ВЕЕ = ALLL ACI NEO EL OI YC Itr ZLC? 
Pirate ПАЛЕ ЕВЕР Ta FING SUAFALE En. 
ПО 949 [=Z 9 
po 15-1. 
ШІ TZXPECI2 
ШЕ Sos vPE I» 
GX1,45zXP4 | FFA 
nc n A IS 
Gel, Вур | ЖҮЕСІ) 


DONC añ a 


3 Есіз-?РСІ) 


DO 953 1=1,6 

GCT, 1226C98,12 

ES I=L, +66, 1) 
GOB, 1 =G07, | +68, 12 
Fo, =F COS 

FCS =F CS FCG) 

Есе = ЕФЕ (Е) 


ac 





Е SOLVE MATRIX EQ. G*CZ=F FOR CZ 
CALL MATINS¢G,9,6,F,6, 1,0ETERN, 10M, INDEX) 
I OC I D EG 12 GG TO (953,9600 IT 
ТЕБЕ - 1 
WRITECO, 954) 

994 FORMAT (33H ERROR IN INFUT - SINGULAR MATRIX > 
LINE=LINE+1 
102=102+1 
50 Т0 450 

955 IT=2 

FINO NEN NORMAL VECTOR 
XN-XN-CZC2 2 *X3-C2C3 )*X4 
ЕЕЕ үз=Г? ез ә*үа 
Z2H=2?H-C 2.2 ES LL CS FE 
SUIZA YH 282 
AN=XN*A 
VN=VN*+A 
ZN=ZN*A 
GG TQ 945 

C STORE DATA 

Gea B¢d+ 1 =xXPC 1d 
с а= чре |> 
BC S43 2=XP C2) 

BC U+4 =F C2 
BC J+5 .=XP C3) 
BOJA EXPO 2 
Еси+т оу ә 
BOE IAS 
ELIAS 
Bv. H. 10 y=23 
BC 11252X4 
E+ 12 =4 
Bt J+ 132-74 
Бо 4140-2610 
EXA SO 
Bé 16852C2€C5) 
BECAME 1? 202662 
IF «K .LT. ?ҰНВУНПЕ) GO TO GES 
10?=107+ 1 
K 1=K+K 1 
K=0 

aes CONTIHEE 

TiK+1 =X 

TOK +2 S40 

TAK4S 9=20 

TCR +4 y=xH 

TES =H 

TCK+6 y=ZH 

TCK+? »-RO 

E COMPUTE QUADRUPOLE MOMENTS 

SAPO IAAF OZ) 

9 РЕЭЛМФВБЫ<5) 

Pere LP CS I4XP C2 > 

AI 4=SXP CS MXP C4 > 

MIS=XP C2 XP C4) 

ALIS PLA AV PCR 

VIZ2=YFL 104 WP úd ) 

ШӘ-УРСЗОЖҰРС2) 

vI4dzVPCIORVPCAD 

VISHYR C2 AYP C4) 

Fi1zH1/24. 


C 





R2=A2/24. 

R3-R0/12. 

ВХХ=СХ1 1**2+Х12*®*д 9*R Le (X | PEZ AX | F2 KAZE HDA 

axvecx prey E 8X L2VI 2 )R eX LV] 2eX LV | 4 AREA | S68 I SERS 

Нүу-сү| жеу | DRAE Le Cy | ekz | door 2 DRE Zr S 
COMPUTE SOLID ANGLE 

XX=XC-XCENTER 

YY=YC-YCENTER 

22=2C-2CENTER 

Хі-ххжха үүжуз%22%23 

YOGA sq ez2*z4 

Z izXX*XHeVVEVNEZZ*ZN 

ШЕ) БИЛ ҒАР 

RCLI=FD#*3 

RSU=BCU**2*FD 

SA=EA+2 IH CADFRCU- CLARKE PERZ+Z DEF EX DER ) 


1 нит татар ле ал кто Т 1S жЕ ТЕМ ЕЕЕ) 
EC. dE SAN. 
BiU (GozRXN 


BH 20052RVN 
ERRIA TESTS 
ASS MECA PCS exe CIA) CVPCSO-VECID2, ü. 
NECE XELA CYPCA4O-VPCZ 02, U. 
FL=. 3 *AMAX1C01, 025 
CEZ3=ABSC СТУ > + HES Bis de 2 
Pare eects treaties) ОТ FLY 001 » GO TO 970 
Ba Ta ІТ. FE*.3 » DO T0 977 
BRIENRITEIG , 0753 С223 
ENSSFOHHHTCZOH QUESTIONABLE FOINT -FOOR FIT ,СЕ14.3) 
| ERRZ 1 
LIME=L 1 ME+ 1 
E EC YFCZ 2 FAPL1dA БЕ. 0. .RHD. 
1 VELCZIEXPLZI-VPISIKXFIZD GE. O. .HBHD. 
2 Есе PAIS POR AROS) .GE. 0. AND. 
3 BER REXE AFE LFF DE. G.) bl TO God 
US HRITESE, 10003 СЕО НЕ O 
1000) FORNAT<41H ERROR IH INFUT - CROSSED 0R CONHCANE OUSO 


1 Aa SA 
IERF= 1 
LIHE=LINE+1 
IDS=1D3+1 
РЕНТЕГЕЕСПОЧСАР(2)-АРСІ»), СҮРСә>-ҮРСІ»2 ПОЖХРСЗ)-ХР(22% 5Ш2ССІРСІ12- 
1 j ri np dy G FH ULI TITIES, YEA YES. 2 


Е ТАШ GT. ЕЕСЕЕЯ2 > GO TO 256 
LINESLINE+ 1 
HRITECH, 901) 


G61 FORMAT (24H WARNING LONG THIN QUAD. > 

EU < 21 .GE. 0. 2 GO TO 35 

eee WRITE (6,55; 

За FORMAT (23H QUESTIONABLE POINT = INWARD NORMAL > 
[ERR 1 
LINE=LINE+ 1 


J. EDITE QUAD INFORMATION 
351 1F CIEDIT1.E0.2.AND. 1ERF.EO.0) 60 TO 354 
IF € HEDIT? .E0. 1 > GO TO 354 
GO TO (256,357,258,359) НО 
356 WRITECE, 519 M, XCF 19, XCP2>, XCPS9, XCP49 XC YN P 6264) 
I М, РТУ Меру Мерз, МСР) Ус ЧН РЫ CZCS) 
2 Р, 2Р1), 2Р2), СРЗ), P4, 2C, ZH, CZC D, CZCEO 
G0 TO 360 
? MRITEC5,510. H,XCF22, XXPZO, X€P40, XCP 2, XC, XN,Rü — C244) 


M C 





Гг 


1 
2 


250 
1 
2 
354 
1 
2 
360 


SI FORMAT CIH ,13 ЕЕІ4 -5/1Х,13,8Е14.5/1Х,15,6Е М4.5/2 


г, 
С 
ya 


$00 
402 
404 


45 


410 


420 


N,VCP20, VCP320, VCF4O, VCP 12, VC, V 


SED 


P,2(P25,Z(P32,2CP4 3, ZCP 12, 2C, ZNL CZ C10, C266) 


GO TO 360 


HRITEC&, 51». H,XCP42, XCP 10, XCPZO, XCP32, XC, XN, Rü 
NYLEG VCP 12, YCPZO, VCESD, VC, VN.FL 


C264) 


ә и Е < 


EID a p S (Ps, 2C Zh СЕС С2СБ 


60 Т0 360 


ЫВ!ТЕ(6,512 M,XCPS), XCP4), XCP 19, XCP22, XC, XN, AD 
N, CES), VCP4 2, VCP 12, VCPZ2, VC, VN. FL. 


C264) 
‚С2(5) 


Р,2<Р32,2(Р42,2(Р12,24Р22,26,28,СР<12,02(6) 


CONTINUE 


LINE=LINE+4 
IF (LINE.LT.50) GO To 354 
WRITE (6,42) 

LINE=O 


4 CONTINLIE 


Ш=. +20 

| =F 
CMC I S=UNF 
F=F+1 
NÜ=ND+1 
fERR=O 


a K=K+? 


К. WRITE GUT BLOCK ОЕ E ARRAY IF FULL 
[F CLT. 240) GO To 400 


S WRITE (042 0,58(12,1-1,2400 


üzP 
=0 
O E TE ULE CUER FT. ARRAY 
CONT | HIE 
niznit2 
Ni=Hi+z 
ШЕ СЕТ СПБ НЕХТ БЕСТІПН 
NS==H5 
09 405 М=ММІМ, ИМАХ 
DO 403 N=HMIN, NMAX 
IEH, N o= 
FÜ=1 
ЕЕ ЕП ОО; 60 TO 410 
МИНХ=Н! 
МНІНеН! 
HMIN=HI 
NMRA | 
GO TO 420 
МИН = | 
MMIN=MI 
HH NSHI 
ЧИЕН | 
NE=ME 
UHRZUN 
НЕГЕ Q) GO TO 50G 
WEITECO6, 452 ICHTREL , NS 
LINE=Ü 
GO TO 300 
HRITE (042 Q,ECID, 121,240) 
NP=CK+K 19777 
ISHP = ISH + 1 
EV IU COS, SOO, SEO, SOD, SHP 
SH=SA+SA 
SA=SA+SA 





590 SA=SA+SA 
01 WRITE PARAMETERS AND T ARRAY ON TAPE 
999 J = 1 
IF CITAPE.EO0.1 ANO, EOFSA .NE.O 3 GU Ta 401 
ogy WS¢d> = X] 


292 
600 
601 


BOS 


519 
ei 
Al 


mee 


zz 


ho po 


cn cn Ch om Cn em Cn 
3 C0 69 62 bh ы 
لا ار‎ — Cy O CG CH 


o 


634 


p 
= Í 
ЗЕ 


ел cn 


24 


СТ 


640 


HS(J+20) = vi 
WS¢J+40 = Z| 

= +} 
ІРСХІжж2 e VpSEZ + 22) 599,599 ,598 
WRITECH,600> ХІ,УІ,21 
FORMAT 11HCEXTRA FLOW, 10X,3F 12.5) 
RERD(S 40» XI, VI, ZI 
IF «EDFCS? EQ. Q GQ TO 597 
= 
YI=0. 
21-0. 
50 TO 597 
IF CISP.LT.O) GO TO 805 
WRITE (03) (FROECID, 121, 150 
WRITE (03) €HSCI), 12 1,2200, IEDIT2, LEDITÀ 
WRITE (OG) <TC), 121,82 
WRITE (032 CDMC13, 1=1,NP3 
CONTINUE 

H1 CHECK SOLID ANGLE 
HRITE (5,50) SH 
FORMAT<14HOSOLIO ANGLE = ,F8.3) 
REHIND 04 
REMINO 03 
103=101+102+103+104+105+1D8+107 
IF C105 60.0 AND. NF.EU.NOüE) GO TO 638 
WAITEC6, 625) 
IF C1D1 .GT.0) WRITE<6,628) 101 
IF (102 .GT.0) WRITECE,629) 102 
IF «ІЗ 07.09 WRITECE, 6209 102 
IF (IO4 .GT. O) WRITEX6,631) 104 
IF ¿105 .GT. 05» HRITEC&, 6325 105 
IF CIDA GT. 0» MRITECB, 6225 IDE 
IF CIO? GT. 03 NRITECR, 834) НЕ, МАХНОЕ 
IF CHF ONE NOE) WRITECE 637) МЕ, НОЕ 
STOP 
FÜRMRTCSSHOÜFRTRL ERRÜR IN DATA - PROGRAM ABORTED? 
FORMAT<1HO, 15,218 QUADRILATERALS NITH ZERO ARER 
FORMATS 1HO,15,44H CUADRILATERALS GENERATE A SINGULAR MATRIX ) 
FORHRTC1HO,15,26H  CROSSED QURDRILRTERRLS ) 
FORMATCIHO,15,32H SECTIONS HAVE TOO MANY POINTS > 
FORMATCING 15,24H POINTS HAVE INUALIO M,N INDICES ) 


з FORMATCIHO,15,52H SECTIONS GO НОТ HAVE CURDS ARRANGED IH GROUPS 0 
Ir 4 ) 


FORMATS 1HO,15,46H QUADRILATERALS GIVEN, EXCEEDING THE LIMIT OF 
14) 
FURMHTCIHO, I5, 28H — QURORILATERALS GIUEN, HOT ‚1427 


ИШТ CIS? .LE.0 со то 540 


WAITECE 6309 ISF 
FÜRMATCTHO ISF= ,14, 19H - PROGRAM ABORTED 2 
STOP 
CONTINUE 
02 RERO PEFS? AND TRANSFER TO IT 
STOP 1 
ENO 


ЕМНЕТТОН 6502(Х5,%,22 





с) 


CCOO (Uy Ú) ú CRESPO O O COS СТ IIND 


myc 


Су сз Су 


COMFUTE SOUAR ROOT OF В*®*®2 


той 


60 
2û 


es 


100 
103 


R= ABSCXO+RBS CY)  +ABSIZO +. DONDDUDODOOOD 1 
IE EZ -2жж? 

R=R+ES /A 

R=.ZI+R+ES ¿E 

R=R+AS / 

SOZ7= .25*R+AS /R 

RETURN 

ЕНГ 


SUBROUTINE MATINSCR,NR,NT,E,NC,M1,DETERM, IC, INDEX > 
MATRIX INVERSION HITH ACCOMPANYING SOLUTION OF SIMUL. El. 
PIVOT METHOD 
FORTRAN IU SINGLE PRECISION NITH ADJUSTAELE DIMENSION 
NOVEMBER 1971 5 2000 NAVAL SHIE A & O CENTER 
WHERE CALLING PROGRAM MUST INCLUDE 
DIMENSION ACNE HAO, BONR NC, INDEXCNR, 3 
WHERE HR, NE ARE DIMENSIONS OF A,B, INDEX 
Ni IS THE ORDER OF A 
M1 IS THE NUMBER OF COLUMN VECTORS IN E МАЧ BE 0) 
DETERM HILL CONTAIN DETERMINANT ON EXIT 
I WILL BE SET BY ROUTINE TO 2 1F MATRIX A IS 
SINGULAR, 1 1F INVERSION HAS SUCCESSFUL 
MATRIX A CINPUT MATRIX? WILL BE REPLACED EY A INU 
MATRIX E: TRE COLUMN VECTORS WILL BE REPLACED 
РУ CORRESPONDING SOLUTION VECTORS 
INDEX: NORKING STORAGE ARRAY 
IF IT IS DESIREO TO SCALE, TRE DETERMINANT CARO 20 MAV BE 
DELETED AND ПЕТЕР РРЕЗЕТ BEFORE ENTERING THE ROUTINE 


ПЕНЕН АСМА. НЕЬ EINE NCS, INDEXCNR,3> 
EQUIVALENCE («ІБОН,«ҒОЫ), SICOLUM, JCOLUM>, CAMAX, T, SHAP > 


SRL ZATIGN 


МеН1 
И=ИЇ 

DETERM=1.0 
Dü 20 J=1,N 


INDEXCJ, 3220 
DO 550 [z21,H 


Sembee FUR PIVOT ELEMENT 


AAR = 0.0 

Co 105 J=1,4 

ШЕІЗШЕЛІІ! 22-1) GO, 103, 60 

00 100 K=1,4 

PME NCK 29-19 30, 100, 715 

Ег AIRI -ABS CAC Rod 8S, 100, 100 
{POW = J 

ІСОМ = K 

AMAX = RES CACI Koo 

CONT | HUE 

CONT {HUE 

IHUEXCICOLUM, 32 z iurc DEP Soi 
IHDEX£1,15 = ТРОН 

INDEXCI,25 - ICOLUM 


INTERCHANGE ROWS TO FUT PIVOT ELEMENT ОМ DIAGONAL 


Сі 
Cal 





су Сус) 


mac 


C) 633 €? 


140 


400 


13 


ánÜ 
RR 


ot 


-J 
= 0 
ІҢ 


210 
720 
730 
& 10 
715 


IF «IROH-ICOLUM? 140, 310, 140 
DETERM- -DETERM 

DO 200 L=1,N 

SWAPS ALIRON, LO 
ACIRON, LO=ACICOLUM, Lo 
ACI COLUM, L =SWAF 
ІЕМ) 210, 310, 210 
DO 250 L=1,M 
SWAP=BCIROW,L> 
BCIROW, L=BC COLUM, L> 
BCI COLUM, L =SHAF 


DIUIDE FIVOT RON EY PIUDT ELEMENT 


PIVOT = ACIECOLUM, ICOLUMO 
DETERM=DETERMYF | UOT 


I ACICOLUM, EOLUS = 1.0 


OO 390 Lz1,H 
R£ICOLUM,L X RCICOLUM, L2 /F LUOT 
ІЕ (М9 350, 350, 360 

ШІ 500 L=1,M 

Bé COLLUM, L)SBCICOLUM, L УРІМӘТ 


REDUCE NüH-PIUDOT ROS 


ШЕН SSC LIST, 


um 


IFCLI-ICOLUM 400, 550, 400 
T=AtL1, ICOLUM> 

ACL1, COLUM =0.0 

DQ 450 1-1,М 
ACL1,LOSACL 1, LO-AC I COLUM, LT 
IF (М) 550, 550, 460 

00 500 1-1,М 

B4L 1, L2SEXL 1, L-BCICOLUM, L 2*T 
CONT | HUE 


с. 


ІМТЕЯСНАНСЕ COLUMNS 


00 710 1=1,N 
L=N+ 1-1 
IF CiNDEXL, 19-INDEXCL, 239 820, 710, 


ott 
tat 


 JRON= INDEXXL, 1) 


JCOLLIM=INDEYXCL,2) 

ПО 705 К-1,М 

SWAP = ACK, JOU: 

ACK, JROW SACK, JCOLUMD 
АСК, .ICOLLI 28 SURP 


CONT TRUE 


CONTINUE 

UO Tal Ks, 

Ж ЕСІНВЕЛСЕ, 32-12 715, 720, 715 
CONT I NUE 

CONTINUE 

I2 1 

БЕТПЕН 

ІБ-2 

50 TO 310 

EHO 





APPENDIX II - XYZPE SECTION РЕ2 


PROGRAM PE2COUTPUT= 128, TAFEG=OUTPUT, TAPEOS, TAPEOZ, TAPEOS, 
TAPEOO, TAFPEO1,TAPE11,TAPEOS, TAPE4=TAFEOS, 
TAPE1=TAFEO 1, 

TAPES=TAFEOS, TAPE2=TAPEOZ, TAPES=TAPEOS, TAPES=TAPEDO) 


hà 60 — 


XYZ POTENTIAL FLON PROGRAM.  UERSION 4 SECTION 2 
БЕНЕШТЕС NATRIX COEFFICIENTS 


c C C) C) 


COMMON Ex291)3, 198000, Y 1010003, USC 10009, VUSCIODOS, EI auro, Czc 9 
200), C30 шой), HSC300),PROBC1S), 
SULLE, VYS), UZC) 

ІШІП! ЕМГЕСУЗ Ү2) ПОЕ SC O ЕМІП, КЕСТІ» ҚМО 
ОО 52125 РЕ) _‚‹ы5‹20@), 1Р5) ‚сЕМ,МКУ 

INTEGER SN 

БІК-1.0 

НЕЕ 

БЕН 032 CPRUBC 13 1=1, 153 

URITECE,S) 

FüRMHTCAaHOXVYZ POTENTIAL FLOW PROGRAM SECTION 2, VERSION 4 2 

FORHATCIHU, 1594) 

WRITE <6 ,90) (PROBCLO, l=1, 15) 

Е A. EAD PARANETERS, T ARRAY, FIRST BLOCK OF E ARRAY 
RERDCOS > MSc PX I 
RERD(OS3 — CTCIES F8 ME) 
REACICOS > (ECT, Iz1,241? 
Ë. START COUP OUER GURÜRRLRTERRLS 
Қ-1 
ШЕ | 
Р-1 
К = | 
USE 
NPH=S+¢ CNP +4 oS 
КИМ ЕНЕР 1 

ІЕСЕС 12-Р2505,205,595 
FORMAT POINTS QUT OF ORDER  Ec17=",1F4.0," F=",1F4,0) 
dez 
XIB) 

МісЕ4М%12 
Х2=Е 1+2 3 
essi: 
X3=E (+4) 
C Waste 
X4=ECU+5 9 
Y4=E +Û 7 
XN=TCK+4+3 > 
VHET CK +4 0 
ZII=TCK+5) 
XCZTCOE 
YCETCOR 12 
2С=Т«К+2› 
A=T {E+E ? 
AX-TBGHM?) 
A= +R 
EX=BCI+9 9 
XVzBG 102) 
YY=E ¢+ 1 1 
Zr=Bil+12) 
С Ci COMFUTE LENGTH OF SIDES OF QUAD 
Di2-SQzFCX1, X2, V 1, 2,0. ,0.) 
D23-S02F (X2, X3, V2, V2, .0, .0) 


uo 
Cg 


у 


hi 
Ci am od 


Пі ың (тал 


Ra pa 


LO 
Cn 





34-502 (X3, X4 , V3, V4, .0, .05 

D4 1=S02F X4, X1, V4, V1, .D, .0D 
C2 COMPUTE SLOPE OF SIDES 

ГЕСХ2-Х302505,300,505 

CIES. 

60 70 310 

AOE AN 23D 

CI23-0. 

IF(X3-X42315,311,215 

С134= 1. 

Go TO 320 

СН24-(Ұ4-У22/0Х4-Х2) 

С1з4=й. 


Ш ІГСА4-512522,521,222 


СІ4 ті. 
Bu ToO ss 
CH4 12 CV 1-45 / CX 1- X4 » 
MESLIN 
IFCX1-X2 7335, 331, 335 
"zl. 
бо Та 340 

AAA A PS 
СІ 12-0. 

C3 COMPUTE QURDERFOLE MOMENTS 
CIXX-ECUOT 172 
CIXVZBG JA 182 
CIYY=E J+ 197 
CV 1220.0 
CX 1220. 
CY SS=Li . 
Cx23=0. 
Crad4=0. 
CXS4=0 . 
ЕРЕ! n 
СХ41-0.0 

ШІ COMPUTE SIN AND OG ГЕ eee reeset ERLE Silu 
ІР (П122 8341,4242,9241 
СУ12= (2-9 19/6112 
EX 1Z=CX1-XZ 37012 
IFéDz2 28242,0244 9545 


Cc) ccr 


343 емуезтіұ2- V2) /02 3 


CX23-CC2-7 X3 2MDz3 
І(ПЗ428345,9348,9245 
Cy24=(v4-vS/D34 
CX34- (X3-X4 2/D34 
IF(D4 120347, 9348, 0247 
См 12 CV 1-45 /D4 1 
CX4 12 OX4-X 12/04 1 
C5 COMPUTE NAX LENGTH OF QUAD 


у 


249 5Т=НЕ КЖ 1-ЖАЗ 


ST2-502F Q2, X4, V2, 4, .0, .0) 
STEANAX 1(ST, STZ, 012, 02: 370346041) 
D. START LOOP OVER NULL PTS. 


(f) Cf) 


43 1=K0 


MILLE 

L=1 

ІЕСІР5) 95Б0,8360,985350 
I (IPS) 4355, 93512, 9392 
EC IFF 9360, 9360, 9535 


ӘЗ ӘСЕ 
Caf JC =. C 
C3CAC »z. 0 


О 
Ti 





сә С 


ХС0=Т<1) 
VCOSTCIE 1) 
2Сй=Т1+2›) 
ХНй=Т(1+3:2 
YHCET C144) 
ZNO=TCI+5) 

E. COMPUTE DISTANCE BETWEEN QUAD AND NULL FT. 

_ DETERMIN METHOD 

RPO=SOZFOXO, XCO, YO, VCO, 20,200) 
IFCRPO-ST#4 3250, 250, 460 
Х=(ХС0-ХСоЖХХ+(УСП-МС ож УХ+(2СО-2Со*2Х 
VzCXCü-XC XXXV CYCO- VC Y v cZCü-zZC 2424 
Z- XCU-XC YXH CI CO- VC YvHe ZH 
IFCRFü-ST*2.03355,355,400 

F. COMPUTE VELOCITY COEF. BV EXACT METHOD 


255 RieSü2FCX, X1, V, V1, 2, 0.) 


tat 
Ті 


2 Ë 3 


Eu asHIRHCCCHZS*E 


BELEK WV v2 2 cl. 

ES HORDE XS V ovs s. 
H42S02FCX,X4, V, v4, 2,0. 
IF ««Бізгге ДЕ, Biz» Go TO 1000 

IF C£CRZ+RS)> .LE. 0235 GO TG 1000 

IF CCRS+A4) ДЕ. Пз4› GO TO 1000 

IF CCR4+h1> LE. D41) GO To 1000 

CLRisRLÜüGZCR T4R2-D 12 ^ / £P T4R240112 22 
CLRZ-RHLOGGCCRZ*R3-DZ3 3/(IZER34U0Z3 22 
CLAS=ALOGC CRS+A4-D24 5 / CAS4+AS+054 9D 

CLAF=ALOG* CR4+F 1-04 fo CR4+h 1+04 133 

TUXSCY 12 *CLR 14 CV ZS *CLRHZ«CV4 "CLR CVAÀ ТЕСЕН 

TUYZCX12*CLR 14 CXZzFCLBZ 4 CX S4 *CLRS4 0X4 PFCLRA 

TUZ=O, 

TE aS 010 375,361,361 

zS0=2+#2 

E 122204 CX-X 1 62 

ES ¿Sue X= e EZ 
O E ЖЕР 

Eg= 2504 O ER 

H1= CY AX 

MES 2 YE CX-X2 > 

H3=CY-YS F X-X3 ) 

На4-(Ұ-уүджах-хао 

ІРІСІ 122363, 363,354 

Hg 1zCCH12*E 1-H 15 /CZ*E 12 

EEC-OLCMIZUTEZ-HZOZCZ*RBZO 

RT 12RTRRCHE 12 

AHTZ=ATAHCHSZ 5 

TUZ=RT1-RT2 

[FeCl 23 9366, 366,367 

2-H2 yC Z*R2 12 

RT4sRTRHCCCHZS*ES-HZO/CZTRE32) 
TUZ-TUZ«*RT-RTAÀ 


"© "` "s 


ШИГЕСІЗ4 2565,306, 380 


RTS=RTRHCCCH34*E2-H2)/C22#R2)) 
RTESRTRNCCCH34*E4-H4 )/CZ*F4 2) 
TVZSTUZ+ATS-ATE 
1F(C141)370,370,375 
RT?SRTRHCCCHA T*E4-H4 2 CZ4F4 2) 
RT&-RTRHCOMA TAE -H 3 2F 199 
TUZ=TUZ+RT?-ATS 
90 TO 450 
G. COMPUTE VELOCITY COEF. EY GUADRAROLE METHOD 


A 
==] 





400 RFOS=REO**S 


461 


470 


MT 
gan 


451 


ВРО?=СВРПЗ**2 У*ВРО 
HS1= Х/НРОЗ 

XSU= EL 

YEY 

2<0-2 жж? 

ip or XX 

05-Х50%250-4.жҮ50 

|gzzxXx*09  *PS«30 . *xS0»/RFO? 

|323 "ipe /RP? 

HS4z3 .x»rüs /RPG? 

TUXZR*US 1-C | X NELIES-C IL XX*USZ- CIN SELIEA 
Hg 12V /FPLES 

WS2-V*ca O82 YEO) RFA? 

TEENS I ZEN ENSS-EI RT ENSI- CINE 


Бет REES IE IAEA EA EDS PPO) 


UA ISI TILDA TUVE TUZ ЖАН 
UY |Sy=TUX*VxX+TUV*yvy+TUZ+*VH 
E TUE AA TUVE TULF 
GO TO 470 

Berle VELOL ITTY COEF. БУ MONOFOLE METHOD 
HFFIISZH/CEPLUES S 0 
UKE SS=  (OXCU-XC 2t nni Fs 
мус = QVEQ- VC PÓREFOS 
UziiSo2  CECU-¿CIFHAPOS 

O TELET ET PT FH PLANE ПРЕ SVHETRY 
GO TO¢4s0, 455, 490, 495, 300, B05 AI 51585 I“ 

DIRE SERIE OG EOBHEESUEE Gr TROUE REGISTERS 

dizdi 
«т 
ONAN (19 
UDY=UYC1) 
UDZ=U2 1) 
NRT SUNET 17 
J2J1 9=UX( 15 
Uae 14 3=UXúC 13 
Nitra) 
UZEI pay 17 
Чох 1+ 0с 
U1CUTe222U2615 
U20J41+2)=UZ/ 1) 
acl {+2 ші? С 15 
IFCSYM) 530,530,481 


IS=2 

Ке SYHEIFY 
Y= 
GG TG 345, 


IFCEUM-19517,517,486 
ХУ SUMETRY 

Е 

200=-2 

Ga TO 345 

|54 


у= Сй 
GO TO 345 


E ii2 310,310, 42% 
VZ SYHETRY 

|5=5 

БЕЙ--ХС0 

Ga Ta 245 

Ier 


“Cd=-' co 





сас 


515 


315 


cn 
y 


276 


зей 


GO TO 243 y 
{5=? 
éCQ=-200 
GO Ta 345 
|Sz8 
vVCü-- VC 
GO TO 345 
J. ADD CONTRIBUTIONS OF ALL REFLECTIONS 
U1CJ122U 1C 104UX COO*UX C? Od UX CB Ux cS 2 
U26152U2C.12-UX CB MIX C? MKUX CS 2- Ux C5 2 
Vac IS URAC IA A AA E ORUXCS 2 
ТІНІ 101% 12-МУ ССС НА АБ 2-МҮ С» 
MECA т у= уг + нус IH DIVE NT) 
ERSTE DENE 2-UYC? FUN CR -UY CI J 
O A A AO LA I ) 
AID IA IDA AR WEIS} 
ITE SUSIE ZES IHE y+U2C6 +U CS) 
UEILS IR UA 2 2 
ҚУАТТАУ РЫ <р асы 
USCUI152Us3CU12-UXCd -UXC2)5 
Micdietosuté Ie 120«Uy eau c3) 
Uzct die Tosuzcdd ie 10*UVe4 06M cz p 
OCI 1+ 1 ӨРІМ ЫЗ 
Ulcdiezosutcuditezo-Uzcds-Uzcs) 
VCS 142 UEC 1625-U2CA M UZ CS) 
USCJiv25-U3CUI T4242 HILO FULTS 
UIC dt SU Cd Os X Cz 
U2CJ15zU26125-UXC25 
TSCA EISI TAIX ә 
UTCdTeiosuicudie15-UVCZ2) 
VZJ {+ SNZ 016 10€*UV C25 
UoCd Tt 15o-UsCd Te 12707022 
PCI EZ =I E IE ZO AZ ZO 
AI A IA ) 
UscudiezosuscdiezMUuzicz) 
СС =XNOFU 1 CUL 10 VRNQURU 16.1 18 12€ 2HCPRU 1 0.1 192 5 
DEzdz)-XHOvU2c. 15-2 C1 15 2Z HO u2 6.11922 
C362 )—XHOWUS CU 1 e VROPUSC EUR T9 ZNOIFUG C E82 2 


y Ms US 


JC = IC + 1 
TASTE CUESFICIEHTS 
KR. WRITE COEF. ON THRE OP DRUM IF STORASE AREA IS FULL 
geese (0013570, 555,555 
a= 
VIK ELK 
Uz6 1328L&K 
Мз MOBLE 
IFCELK-626.0) 560,563 566 


| ИВ1ТЕ (01) ELK,U1,U2 US 


GI TO 562 
КЕНІН 01 
ARITEC11) ELK,U1,U2 3 
ELK=ELK+1. 
IFC JC-Sti 1 550,37 1,571 
1OW=1 D+ 1 
WRITE «O75 IDU,CI 
WEITE £025 IDII, CZ 
WRITE (09) 104,02 
«С-1 
KO=KQ+? 
L=L+ 1 

L. ENG OF LOOF OVER NULL PTS. 


99 





C»? C) 


58 


— 


620 
SFA 


Aal 


1000 


MW 
«i егі г ! Ü 


IF£KÜ-KM 2343,345,55 
Сіс =0 
C2CC 520 
C3c.C 520 
ІЕСКО-КММ 0541, 565,585 
Р-Р4 1 
K=K+7 
420 
IFCK-KM 596,586,600 

N. END OF LOOP OVER QUADS. 

H. RERD NEXT BLOCK OF E ARRAY IF NEEDEOD 
IFCJ-24 13206, 520, 500 
RERDCDUA XCECIT3, 11, 24 02 
= 
ІЕ<ЕСІЗ-Р2505,2 
WRITE (6,98) Ec 
STOF 
IFCBLK-628.0) 610,620,630 

0. MAITE REMAINING COEF. GN TAFE OR DRUM 
URITECO1ELK,U1,U2, US 
REMIND 01 
DÚ Tü اجج‎ 
REWIND C} 
MRITECTI? ELK,UT,UZ, US 
REWIND 11 
WRITE (92) IDH, C 
ПЕТЕ COS? | HUN, C2 
HAITE (02) IODH, CS 
БЕМІНІ UZ 
REWIND 03 
БЕМІМІ C4 
БЕМІНІ 0S 
FEUIND 0% 

ШЕ ARS FER TO PERS. 
Go Ta Soon 
ЕШ>ШТЕСЕ zu00) L P 
ШЕЕНТКЗН І- 15 205 Н F= ,FS. 12 
CONTINUE ЕГІ | 
ФТОР 2 
END 
ШЕТІН SUZFCA 1, 42,V1,V2,21,22) 
=; 


ip | 


u = 


nt 


2-21-22 
BN e 
HB-RERCXORRHESRCVORHER(QZÓR 1.0E-20 
F=F+Fs zF 
ES. ZIPREAS ZR 
F= Е+Е5 (5 


t nios 
SüzF- .Z5S*BTRS /P 
RETLIAN 
END 
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APPENDIX lili — XVZPE SECTION PF3 


2 TAPES=TAPEOS , TAFEG=TAPECG , TAPFES=TAFEOS 3 
С 
C XYZ POTENTIAL FLOW PROGRAM VERSION 1 SECTION 3 
С SOLVES MATRIX EQUATION FOR SOURCE DENSITY 
E 
COMMON 5М‹654 >, ИІР 2,55, ВЕЙ >, РАСОВ 15 >, 2220 >, 0МС650 >, 
1 BC2205», COEFCIOOD, XNCESOD, YNIRSOD, ZUEGSÖD 
EQUIVALENCE CHS¢21S> EPS) ‚CKK,BCZONDD 
EQUIVALENCE &MIX,H8C20220, XHIV, H8C20820, CHIZ, HSC207? 5) 
EBDDIURCENCE C(HSC2012, NP 2, CNSC208 9, IPS 3, CHSC2 122, IPF 5, CHE C2112, Kl 2, 
1:КМ МЕ. У 
Ss FORMATC4IHIXVZ POTENTIAL FLOW PROGRAM SECTION 2, VERSION 4 3 
ЫБІТЕ (6,52 
READ Mestre 1, 11,139 
HAITE € Û, TOO 129€PROBCIO, Iz 1, 192 
1001 FORHRT& 1HÀ , 1584» 
READ (032 KHS£15,121,2205, IEDITS, IERI Td 
READ СЕЗСЕ ЕКЕ ЕКЕ, СУ VNECI S, 2ZNCIS, SKIP, I1, NP2 
[i2- 523. 14 1590205 
READ (OOOO o, 1=1, NP) 
К. 1= 1 
K2=NP 
240 FORNAT C1&HOCHANGES IN PROB -, 15842 
Be 32270 1220. 1221 
221 RES К т ЕК УК r IS) 
WRITE (6,240 KBiE 2, Ks1, 190 
БЕНО 1225 ЕК5, КЕ 1, 220) 
READ C12 SKIP 
READ CIZOSKIP 
КІзіІР5 
K2=1FF 
8 R. SET CONGITIONS FÜR FLOW üF -1 IH X DIRECTION 
1220 FX=-1 
FY=0 
Fz-ü 
HF-- 1 
С Б. COMPUTE INITIAL REFEOXIMATION TO THE SOURCE 
1240 DO 1250 K=1,NP 
MIPOROZXHCK PFFXEVNCR YFFSRZHNHCOK )FFZ-DNCE 2 
sco KOoz-UIPZE »*. 11326 
E ШӘССЕТ ГЕНТІНІГ SUM VECTOR ТО ZERG 
1280 SHER =U. 
SNCNF+ 1 е0. 
SNCNP+2 2=0. 
SHCNP+3 9=0 . 
SHCNP+4 9=0 
HRITECD, 9497) FX Еу Е 
WRITE CE, Gaza 
9335 FUEMATIZ HON TERAT | OH ES НЫ ІНІ HBI. TUX, THEZJ 
IT=1 
[c=5 
E iPS) 1260, 1260, 1255 
ШЕСІ БЕНОІ 122 С565,К›, K=1,KK) 
DO 1256 K=1,KK 
00 12568 I=1,4 
1256 $61,К›=5$65,К› 
E D. START ITERATION 
1200 EANO=0 


PROGRAM PFP3CQOUTPUTs 128, TRPEO2, TRPEOG, TRPEOG, 
I TAPE 12, TAPEOS, TAPEG=0UTFUT, TAPE2=TAFEOZ, 


101 





1264 
E 
C 
[E 
Ë 

im 

67 

бо 

Ба 

(ti 
C 

S0 
E 
с 

1290 


160 


IF (NE 21261, 1262, 1263 
READ 40221DH, COEF 
GQ TO 1264 


2 READ (089104, COEF 


00 ТО 1264 
READ ¢09>10W, COEF 
ШЕШ 
D. READ FIRST BLOCK OF COEF 
E. START LOOP OUER QUADS. 
DO 1290 K=1,NP 
F. PICK UF SOURCE DENSITY 
SF=SCIC,K) 
G. START LOOF OVER NULL FTS. 
00 1297 КР=1,МЕ,5 
IFC J-G00 80, 65, 65 


ІР CHF G7, 65,59 


READ KOZ) DH, CHEF 


GO TO 70 
|: RERO. COS 3104, COEF 
ой т0 70 
PEAD (GOJON CELS 
«it 
Н. COMPUTE PARTIAL SUM FOR NEXT S PTS. 
SNCKP J=SNCKP O+COEF CJ+ 1 ЖЕР 


SHCKP+ 1 =SNCKP+ 19¢4CGEF CUt+2 OSE 
SHCKP*222SHCEF2 HFLÜEF C3 FSP 
SHE KFS S=SHC RP +S D+ COEF (+4 FSF 
SHCKP +4 D=SNCKP+4 J+4COEF (J+ ESF 
f= 145 
К ЕИН ЕЕЕ UVER MULL ЕТЕ. 
К. ENG GF СОСЕ GUER QUADS. 
ary WAVE 
L. COMPUTE NEM SOURCE 
IF NF? 31,92, 93 
КЕНІНО 02 
GO TO 94 


' REHIND 08 


bÜ TÜ 94 


3 REHIND 03 


Añs== 1.0 

SUHzÜ . 

fi 100 K-ET,RK2 

RIA 3*0 
ШЕ-І-ПЕССЕНКЕ2-ЕСІС,Е>) 
SUM=SUM+TEST 

EE En Dl. EFS) PASS=-1.U0 
CONTINUE 

r ii S ЕШ 1.02 GO TO 180 
ШО СЕ МІХ» GO TO 180 

ШОО ЕШ [з E9. 0> WEITECG, 90) IT, SUM 
IT=1T+1 

ІС-іСс-1 

ЕГІС EC. 0» G^ TO 120 

Cu 110 K=k1,k2 

SCIC, KÓOZENCK D 

SHcK OO. 

GO TO 1250 


|Н. 


B1=0. 
BZ=Ü. 
DA=O. 
01=0. 


— 
= 
Po 





140 


D2-0. 
D) 140 K= K1,K2 
П6а-2ж5(1,К2-5М(К2-5(2,К2 
ІЕСрес GT. Q.) GO TU 122 
Н-Н%5(2,Қ2-551,К2 
DiRZDRH-L'Su 
са Та 125 
A =A %2(1,К2-5С2,К) 
DR-DH«DSG 
D8 12964, K »-S 
ШЕ/--(2,К2-5 
OS3=0S 1- -Ds2? 
[i522 8 6Z K5- SCIT, K) 
DST= 32-058 
DS6=05 1-055 
OS4=DS2-S0 1, KOFENCK OD 
US ?= сзж054- -DSS*DS6 
OSs=05640S5-0544053 
(FeO? .GT. 0.5 GO TŪ 125 
E 1=E 1-DS 405 4+ 032: *D0S0 
D1-zD1-Dz? 
DÜ TO 130 
B126 14+051¥054-D22+40S5f 
01=01+0S7 
CS DES GI. 0.2 GO TO 132 
ЕНЕСІ eo tU DS 
202-055 
GO TU 140 
TERR OS HOSS OEA OS 
o2=02+ 058 
CONT [HUE 
Н-Н/ЛІН 
Е 1=6 1/01 
p2-E2/D2 
ГЕЛ ЕВ. БУ Ga Ta 153 
NC 

HERBS CRED 
IF CAR .GT. .üz» GU TO 145 
00 145 KEK 1 K2 
SCI, K I= =AYCSHCKIRSCL E0480 1 
SHE O=0. 
МЕТТЕКБ, БІН); 


“2,K> 
(2,К2 


EIO FOSMATCZAX, 12HR EXTRAPULAT | 


ЕП ТП ТЕП 


: ББігЕ1-Е21 


EE T-REZCBE 12 

EB 1=50, EE 1 
A A 
BEEZ-REZCEBZ) 

BR2=50 .+BB? 
BBESABSCER 1) + ABSCEP) 


E 


Uh 


IF (5 (ЕРІ .5Т. BEB) OF. <BE? GT. BEE) > GO TO 
ЕШ 150 К=К1,К2 

SCS KISSC2 KE PERCE CT K-S02 Kee BeeCSNCK -S C2, K0) 
SHCK =. 

HEITEC5, 7000? 


FORMATC29X, 17HE EXTRAPOLATION 


GO TO 150 


= DO 158 K=K1, Ke 


SCI KIESNEK?) 


Se Shao-o. 


[C75 


) 


joe 


(mE 


ымы 





ЕС Бі 52 
161 FORMATCZOX,3E12.3) 
AS= 
BS1=B1 
Б52-Б2 
50 Та 1280 
180 HRITEC6, G9» IT, SUM 
Roe K- K2 
182 SC1, KOSSNCKO 
HRITECO3>  €ECL,KD, Ke 1, NP) 
QQ FORHATODX, I3, E 18,82 
007 FORMAT (13HO X UELDOCITVE,F4. 1, 15H — Y VELOCITW=,F3. 1, 
1154 — 2 UELOCITVS,F4, 12 
IFCFZ ACC, 1290, 1400 
E F1 IF THIS WAS NOT LAST FLOW, SET FOR NEXT FLOW 
1390 FZ=FY 
FY=FX 
ED 
MIX=MIY 
MivsHIZ 
NFZHF- 1 
ЕШ ТІ 1240 
1400 REHIND G2 
C P2 RERD IN PFFS4 BHD TRANSFER TO IT 
STOP 3 
END 





APPENDIX IU — X¥ZPF SECTION PF4 


OO OO С 


Nn 


100 


ve 


102 


105 


m 


РЕООВЯЧМ РЕРАСОЦТРИТ, ТВРЕЕ-ОМТРИТ, ТНРЕСОЗ, ТНРЕО1,ТНРЕ11, 
TAPES=TRFEQS , TAFE 1=TRFEU 1 2 


N¥e POTENTIAL FLOW PROGRAM VERSION 1 SELTION 4 
COMPUTES VELOCITIES AND PRESSURE COEFFICIENTS FOR 
Bes Om THE BODY 


COMMON UX 16650) UN IERTD 3, VZ ICHS ua, VCE, CCBAU, E 


om m ө oe = 


RO KIQG 05 USCS, VE SOS EE 0), = Call 550690) 
, ACCU), 6500, 20650). Т4(6502, Т5(650), TELES iQ 
CESS 


DIMENSION’ FRÜBG 153, N5(2203, CU 11000 3, CUECA, CUS CACAO 
EQUIVALENCE WSO RED CWE es КМУ, (HSC2173, Dylan ze Цу, 


HECI, VEL, (MEC208», IPS), (MSC2 12), IFF) 


EQUIVALENCE (MIX, WS¢205), (NIV, WS¢2069, (MIZ, WECZO7 2D 
FORMATC49HOXVZ POTENTIAL FLOW PROGRAM SECTION 4, UERSION 4 > 
WRITE (5,5) 
FEAD (Q2 XPROBCI 5,121, 15) 
FORMATCIM , 1584) 
NF TE € 6, 10D XFRÜECI *, 121, 152 
A. READ PARAMETERS AND SOURCE 
БЕНП (032 CHSC1), 12 1,2202, IEDIT2, IEDITÀ 


READ (0325У612,9С19,2612,745129,Т5612,Т8612,<КІР, 1=1, МР2 
READ CO CONE LD, P= 1 NPD 
П--,5/2. 141520265 
FERD (D33€81C12, 12 1, NPO 
FERD COS MS2C19, 1=1, NED 
READ (O39 SEC, 1=1, NFO 
ШЕН 
Eis 
K2zNP 
ІЕСІРС 108, 108, 12 
Ki=1PS 
K2= IPF 
EBR=1. 

Б. READ FIRST BLOCK OF COEF. 
ІТАРЕ=0 1 


READ (01) ES,CU1,CU2,CUS 
ІЕ «БЕ-БЕБ) 300, 120,300 


2 ÛU 125!=К1,К2 
ШІКІЗ--1.0 s G 2U 
М 1С] )= Seiler Tact} /ü 
Ve {ct a= Soe TEC > ҮП 
Ше л ==) Тас LO 
a E SS рЫ ә /U 
2ےل‎ | = a Тас zn 
UXSC | (= э |)*7Т461) гв 
Іс a= К Aa ¿O 
UZaC | 9=-1.0 Soe pt Gk) eel 
C. SET UF LOOF QUER QUADS. 
«Га? 
BZ EIER UF SOURCE 


51-5 16410 
224552041) 
53552049 
БОО СЕТ UP LOOP DUER NULL PTS. 
OO 160 «Р-к1,К2 
F. COMPUTE PARTIAL SUM FOR 2 COMPONENTS GF 3 VELOCITIES 
UA ICIP OSUI LIPAS 1 Су) ICIC) 
UY TIP SUYA CARS ABRO UCA 


ІШЕ 





e 


Z2 CIP 3SUZ LC UP +S TUFCU TCUC4 2 ) 
UX2 CP )-UX2 CUP )* S2. CU2 C JC N 
UV2 CUP )eUVZ CUP $2.4 *CU2 UC 12 
UZ2(JP )EUZ2CJP )4S2J* CUZCJC+2) 
WS CIP ISUXS CUP SSK CUS CIC 5 
UV SCIP SUYSCUP AS ACUC ORO 
UZSCIP )-UZS CUP S34 CUR CUCA 2 7 
Ces 
G. READ NORE COEF. IF NEEDED. 
IF &JC- 1000180, 135, 135 
35 JC=2 
IFCEER-625.0) 140, 150, 155 
140 БЕНО (01? EE, CUT, CUZ, C3 
GO Tü 160 
150 REHIND (1 
155 RERD (11? EE,CUT,CUZ, CUS 
160 BER=ESR+1. 
IF «EER-BB? 200, 180, 200 
Н. ЕМО OF LOOP OVER NULL PTS. 
180 CONTINUE 
+ل=ل‎ 1 
|. END üF LüüOP DUER OURDS. 
IF(J-NF 2130, 130, 200 
200 IFKEER-B35.05 231,231,232 
231 REWIND O1 
бо TO 233 
232 REMIND 11 
К. EDIT THE VELOCITIES ETC. AND HEITE THEM ON TAFE 


255 WRITE COS IMU UIC VETS ПЕ) 
WRITE C3 2 UX 2€] 5,U24[ 5, U22C1T) 2151, HA, 


WRITE BE CUT TI UEC UZS > mid EO 
229 FORMAT. 1H1, 1S9A4,8H PAGE =, 115) 
REMINC (3 
IPzE 10-43 
| SzK1 
I PRGEz 1 
ПИ ТЕТЕ EU. 1» GO TO 293 
IN dx BE Nu O TO 2865 
424z FüRHMHTCOZHU X FLÜOLL? 
240 ҒОЕМВТСЯН FT., 10X, 1HX, GX, 1HV, OX, HZ, 13X, ZHUX, 8X, 2HUV 8X, ZHUZ , (OX , 
1 SHABS UW 7X,Z2HCP,6X,6HSOURCE, 4X, SHU. NORMRIL ? 
245 FORHRT CIX, I3, 4X, 3F 10, 5, 4X, GF 10.5, 1F 13.5, 2F 11. 5, E 12. 22 
ESUUIFCIF-K2 5253, 222, 260 


П СЕИСОТЕ PRESSURE AND AES. VALUE OF VELOCITY 


255 WRITE (6,235 XFRÜECID, [8 1, 152, IFRGE 
MAITE (6,242) 


HRITE (6,240) 
DO 257 1=1S,1P 
SOS 1€] »P*Z«U0V tC] )*e2 Ug 10] жж? 
UMEC(BECUX TCI D2HRESCUS TCI DORESQUZ CIL 20» 78 
ММ РОМ 
ММ. 258000 ZUM 
ЧИ = SPRCUIMENSO YM 9 
AS 
UNE. — -UXICIOET4CIO — RUVICIONTSCIO — FUZ1C 1 TEE 
ШЕ ШЕПТЕ «6,2455 I, X€13, VCI2,Z€ Io, UX IC DD US ICI D, U2 1€ 1 0, UM 
1 ,CP ,81CD D, UNR 
15=15+50 
|Р=1Р+50 
ІРАСЕ= ІРАБЕ+ 1 
IF(K2-15) 265,260, 250 
260 1Р=К2 


100 





™ 


265 


25? 
270 
Zr 


CORTO 225 


IP=K1+49 

IS=K1 

ENSY СЕС 7) 60 TO 2980 
IFCIP-K2 275, 275, 270 

|P=K2 

WRITE КЕ Zen СЕРПЕ TO 1:1, 159), IFAGE 
WRITE (6,277) 

FORMATCSHO Y FLO 

WRITE (5,240? 

Da 272 I=15,1IF 

USU=UXZ O] ERROR zz zo? 


NM=CHESCHKZO HAES CUYA O) HRESCUEZC IIA U 
UM=UM+USLT ZUM 

Mz. 25 UH 

un = S*CUNFUSO ZUM? 

UNE; ENE FUN E ETS) Јесте 
Er =1.-US50 


278 WRITE (6.245) 1 XC, VOLS ЕСУ UX2CI D, UP2C LS, U22€ LED, Uff 
i CP 52012,04 


те 
202 
225 


£ AL 


231 


[S=1S+ 50 
Р= 1Р0 
ІРНСЕ-ІРНСЕЗ 1 


ІРСІЗ-К22207,207,26П 


IF=K1+49 
І5-Қ 1 
ШЕСІ ЕП» GO 10 Ze: 
IFCIF-KZ 3200. 220,295 
| FZHZ 
WRITE (B, zS3nm x FRÜOECIA, I21, 155, IPHGE 
WEITE (6,2825 
FORMATCSHO £ FLAN? 
HEITE (5,240) 
DO 291 I=1S IP 
VSO=UYX | УЖЖ: яе ружжа ружжо 
UM=CABSCUXSC 0+1 ))4REECUZ St 1 399%, ?O 
LIMzU-US PU 
Unz . 29 *UHe uen lum 
М1 =, DH CUIN+USO ZUM) 
Ер E n 
UNE ou NUI yet CT FMS IPTIDCI?D 
HRITE fe 245) ИАС ООС co DA ON 
I CF ЕСІ» UH 
|= 15+5 
IF=|F+50 
ІЕНСЕ-ІРНСЕЗ 1 
ШЕ(15-Қ22262.2 202 243 
|. CHECK FOR A FOURTH FLOW 
J= | 
Wiel = We 
Uy] 2 Hz: von 
Vel = W&cd+40) 
[FUNC eee peru perz5. 400,400. 301 
IS=K 1 
ІР-К1%40 
EM IL -k25220, 322, 022 
DEDIT Tee VELOCITY AND PRESSURE FOR FOURTH FLOL 
IP=K2 
WEI TE МЕ esa AD OB) |Z! 185 IPRGE 
WRITE Ch SIS UXT UN UZ I 
WRITE (6,3409 
1. COMPUTE VELOCITY AND FRESSURE FOR FOURTH FLOW 


107 





Ln SS i= S IP 
Uxd —OQUODEUX 1€ LUV DRUX?CIOUZ PPUX CI 2 
4 = (QUX DFUV 1C] OUS PRIZE DHL DELIS 0] 5) 
024 -- (UXIFUZ TCIORUN [80226] )RUZ USUZSCI 25 
USO=UX4  **2+UV4 — *o2+4U24 FFL 
UM=CRBSCUX4 J+RESCUNVG ItABS CUZ4 ЕЧ 
UM=UN+USCH UM 
Un=. 25 *UM+USO ¿UN 
UM =. CUM LISO ¿UN 
IF = j.-4USQ5/ UV | F24U2 [FP RU 92) 
333 WAITE — %5<6,3452 !,ХС1І2,ҮС1?,2612,054 ДАУ 0940408 
1 , CP 
[S=1S+50 
IPzIP450 
IFAGE=IFAGE+ 1 
ET ISIC, 222,220 
350 Y = J+! 
60 TO 2904 
S15 FORMATCIGHO ONSET FLOW, UXt=,FG.3, 
S40 FORMATCSH FT. , 10%, HX, AX, HY, Ax, 1HE 
ORARIS N rA REPO 
ile ar Se lu. 3°44 SF U.S, FIS. IFT. SD 
400 CONTINUE 
ба Ta 5000 
300 CONTINUE 
ЫБІТЕ (6,3102 ІТАРЕ, БЕР, ВЕ 
ШШЕСЕННІ СЕНІТНРЕ ,12 17H OUT GF POSSITION,114,6F6. 13 
5000 CONT ! NUE 
STOP 4 
ЕНГІ 


2X, 4HUYI=, FE. 3, 2X, 4HNZI=, FE. 3) 
13X, ZHUX, 8%, 2HUY, 8X, 2HUZ, 10%, 


2 





APPENDIY U_ - XYZPF SECTION РЕ5 


OOOO O 


INT 


C 


PROGRAM PFPFS¿XIHPUT=128, QUTFUT=129 THRPEG2 TRFEG4, 
ITAFES=INFUT, TAPEG=OQUTPUT TRPE2=TRFEQ3,TRFE4=TRFEQ4) 


XYZ POTENTIAL FLOW PROGRAM VERSION 4 SECTION Š 
COMPUTES VELOCITIES AND FRESSURE COEFFICIENTS FOR 
ОҒҒ ЕСПҮ РОІМТЕ 


COMMON  B¢241), XPESCH >, YPESCH I, ZPESOC, WNICSCH), NSC220) 
1, Us 1(5O0Q2, UZ 15005, UXZ C500), Uv/2C5000, U22C5002, UXSCSOQ 5, UV 3 CERO 
2, U22C500), UX C82, UVCB ), U2€82, 8106502, S2C6500, £3 CE500,U 1022, U2 022, US 
3€3) ,PREDECIEO ,XHCBEOD, VNCBSUD, ZNCOBSOO, TREESOO, TX C650? 

4, TV CBS У, TEKEL} 

EQUIVALENCE CKM,WS¢2119), CKM,MK2, CNP WSC20 199, CE HECZ 1020 
1, 672,92), (MIX, WSC20593, (MIN WSC2069>, (HIZ HS(207); 

INTEGER: PAGE 

INTEGER SYN 

S FORMAT<49HOXYZ POTENTIAL FLOWN PROGRAM SECTIONS, VERSION 4 0) 

WRITE (6,5) 

А. RESO INPUT 
RERD(S,252» NÜEF,lEDITS, IRERD 
A. READ THE CFF BODY POINTS 

NOE=HOEF 
DO 10 1=1, НОБ 
READCS 283 ХР), УРАНУ, 2Р1) 

IF СЕСЕ СЕУ „ЕП, 0.) Gü TG 10 
МОЕР-1-1 
WEITECE,G> NOEF, NOE 

3 FüRMRTCTHU,I5,31H OFF BODY POINTS SPECIFIED NOT 130 

ба Та 11 

10 CONTINUE 
11 CONTINUE 


C 


-— 


2o FORMATO S142 


< n 


ШЕ URMATESE 12.5) 


- 


с 
БЕНІ (03) «РЕОВ(12,1-1,159 
WRITE (6,901) CPROBCI3, 121, 159 


S0 FORMATE IHO, 134 > 


HRITECE,81?  MOEP, IEDITS, IRERD 


= р = 


91 FORMATCSHONGEP =, 14 ¿SH IEDITS=,14/5H IREAD =, 14) 


WRITECR, 92) 


FEL IIC pr GFF Boch POINTS / 4H FT., 11%, HX, 12X, 1HY, 12%, 1BZ0 


д 


WRITES, 339 CI, XPClI 2, VPCID,ZPCID, I2 1, NDEP) 


БЕЙШЕЕПІНІС1л,13,2/.3Ғ12.5) 


E. READ TRE PARAMETERS, T ARRAY RHE SOURCE FRÜM TRPE 31 
READ (03)  4HECID, 121,220? 
READ (032 (TXCIO, TVCLD, TZCLO, XHCIO, VNCIO, ZNCIDO, TRCIDO, 1831, НРУ 
READ (03) SKIP 


ШЕЙЕНЕГУ: И5(2202 .£0. 2. 


IESHUESZZUD ‚ЕП. 95. > RERDCOSO SKIF 
READ (033 ¢S101), 121, NFS 
READ (O32 (52013, [2 1, NP? 
READ (OS) (S3(1),l=1,NP) 
С. RERU THE FIRST SLOCK GF THE E ARRAY 
READ (047 (Bx15>,1=1,241) 
К= 1 
Ше | 
DO 100 1=1, NCEP 
D. SET THE FARTIAL VELOCITY TO THE FREE STREAM VELOCITY 





3 


10U 
290 
291 

ae 


Ux Cl o-2- 1.0 
WY ICL =U. 
Uz 1CI5z0. 
UX2C13=0, 
USC) 0=-1.0 
UX3CI =. Ü 
Uv3CI»20.0 
UZ23613=-1.0 
UZ22<1 =0. 
ІЕСВСІ2-Р2291,295,291 
ШЕШЕ Ne SE) EL), PR 
FÜRMATSZEHO POINTS OUT OF ORDER E(15=, 1F4.Q dH P=, 1F4.035 
STOF 
E. START LOOF OVER THE QURDS. 
«-2 


Fi PICK UF QUAD. INFORMATION 


5 ХізБі 2) 


Үі-В ғ 1) 


TZB) 


м2-Е J+) 
Х3-ЕС4%4) 


X4-EC. 32 


e 62 

XCZTXCK 2 

VCETVCK О) 

ZÜ-TZKK? 

AH =TACK OS 

XX-BGCd 2 

YREEC I+ 2 

+ )ا 

ТЕР 

E (+ 1 1 

ZYSBC I+ 129 
F2. COMPUTE LEHGTH OF SIDES OF OURO. 

(iiz-SOZ2FPQX1,X2,V1,V2,0.,0.5 
i22-EGZFCX2, E vs ‚0, w 

Li: ASIERE CAS, Ха, Үл, mu M 1 


F2 COMPUTE ¿Lote E Кіре = 
Az 05, M0, 09 
С123= 
4 TO 310 
A: A 


Ea 
23-0 


ПЕСХа-549215,211,215 
г134-1. 


Gü Та 220 


с С/124-(У94-Ү32/СХ4-Х20 


С i 34=0. 
O E 20 


1 СІЧІтІ. 


GQ TQ 330 
CH4 120v 742 /0x 1-X4? 
Cidi=c. 


0 ІҒ<Х1-Х22335,231,2325 


С112=1. 
GO TO 340 
ЕШ 2-(92-912/Х2-Х12 
СІ12-0. 
РА COMPUTE GQURDRRPOLE MOMENTS 
CIXX=BCJ+ 17) 
СІХУ=ЕСЈ+ 18 2 


110 





92 
+ 


оз 
+. 
r3 


CIYV=B(J+ 19) 
FS COMPUTE SIN AND COS OF SLOPE ANGLE FOR EACH SIDE 
CY 12=¢42-9 19/012 
CY22= (3-2 5/023 
СҮз4-(Ү4-ҮЗ2/З4 
Cv4 12CV 1- V4 2 /D4 1 
CX12=(X1-X2 37012 
CX23=(X2-X3)/023 
CX342 (X3-X4 )/D24 
CX4 12 (X4-X 1 2/D4 1 
Fé COMPUTE NAX DIRGINRL 
ST=SOHFCX1,X2,41,45,0.,0.) 
STZ=SO2F 6X2, 7X4, ¥2, 74, 0, 0 
IFCST-ST2284 1,542, 342 
ST=ST2 
G. START LOOP QUER THE OFF BODY POINTS 
00 530 Jüsi,NDEF 
[Š= 
XCÛ=XP (JÛ 
vCOSVE CU) 
2Сб=2РС 0) 
1-1 


a45 RrFü-SüzFXXC,XCO, VC, vCO, 2C, Z6» 


cc 
„ымыы 


ғ ^» 
CT 
сама 


354 


H. DETERMIM METHOD 
IFCRPU-ET*4 2350, 550, 450 
XscXCO-XC yEXX+ CUCU PY OZONO ZO zx 
vzCXCO-XC Yexv e cv CO- vC vv ca Ca- 2C #2 
Z-CXCU-XC JEANEK e CIO RN CE ССС FENEK 
IFCRFü-ET*z.5»355,355,400 

|. COMPUTE INDUSED VELOCITY BY EXACT METHOD 
В і-502Ғ(Х,Х1,%,71,2,0.2 
RZ=SO2F (X,X2,¥,¥2,2 
AS=SU2F CX X3. V VE, 
R4-SDZFCX, X4, V, v4, Z I 
ІБСӨВІЗЕ23 1Е.П122 60 70 1000 
IFCCKR34R2»5.LE.D22» 60 TO 1000 
IFCCRZ+R4 >. LE. D24) GO TO 1000 
ЕССЕ REO. LE D410 GO TO 1000 
CLR i-RLOGCCR T4R2-D 12 2 /£R £824 D 1222 
CLR2-RLOG(CCRZ4R2-D22 2 /£R24 5240222) 
CLRS-RLOGCCRS4RA-D34 2/CRS4 Rd 4D34 2) 
CLRd-PLDBCCRASR 1-04 157 CR4R 08DA 122 
TUX=CV12ECLA e CVZTSCL Rz e CV S4 CLR CVA TECLRÍ 
TUVSCX 12*CLR Ee CYZSRCLRZACXS4SCLRS«CX4 TSCLR4 
TU2=0. 
IFCRESCZ 3- .ÜC EAST 9375, 88 1, 85.1 


= p w'e 


J 


. 
. 


сэз 
s x 


NI PAF 


4 
J 
2 


y 
af 


. e 


| Z&(ez**2 


ESSAI 

5 E EZ 
ESLE SAS FED 
E422ZS50 CX Xd ree? 

O TD 
AS NA 
SS 

Н4= (0-04 PEGO 

ШЕСЕ| 1223Еб3 303,504 

US 1= CM IEFE R/C 2 FR 10 
WSZ=CCHI2FE2-H2 f ZARZA 
At 1=ATAHCNS 1D 
HTZ-RTRHUZZ 5 

TUZ-HT 1-RT2 
IFCCIZ3»36b,365,307 


111 





C1 


490 


Gu 


456 


490 


АТ3-8ТаН((сМ23ЖЕ2-Н22/(2жҚ2)> 
RT4-RTRNtCCM23*E3-H32/(Z*R32) 
TUZ=TUZ+ATS-AT4 

IFCC134 368, 368,369 


& RITS-RTRNCCCHS4*ES-H2 D / (2832) 


RTB&-RTRNCCCHS4*E4-H4 2 /CZ&FA 2) 
TUZ=TUZ+ATS-ATE 
IFCC1412270,370,375 
RT?-RTRNCCCMA TRE4- HA 2/ CPP D) 
ATR=ATANS (CM4 1#E 1-H 14 Z#*R 19D 
TUZ=TUZ+AT?-ATS 
GO TO 450 

J. COMPUTE INDUSED VELOCITY BY QUADRAPOLE METHOD 
RPüssRPO**3 
RPO?-CRPOS**2 REPO 
NE1= Х/БРОЗ 
XS0=X**2 
Vans y 


PSzVens 280-4 , *xen 
(ssxSso0Ó2250-4 FY SC 
HS22X*(0 *PS4-3D. *XSO/RPFÜ? 
HSasS.*"vY*PS/RFÜT 
NSd=a3 exes “APO? 
TUXESAAHS 1-C I XVENSS—-C INES E—C IVY ALIS 
WS {=Y “RPGS 
ш с=т Ea жүсіэ;ЕРІ? 
TU =R na ii ASS ONES PAZ 
ШЕСЕТТИТЕРЕЗ-с TOC ITI OO трус | уо 2 SAPO? 2 
UNC 1S TURN + TUV ERY + TUSERNCE 9 
IIS STUFF TUYE TULF ERK 9 
SET ТОТ ҮТТЕ ЖЕНЕ?) 
GU TO 40 
K. COMPUTE IHDUSED VELOCITY EY MOMOPOLE METHOD 
НЕРПЗ-Н /ТЕРПЖЖЗ) 
Meets = (ХСП-ХІЗЯНЕРОЭ 
MYC ISo0= ¢CYCO-VC >*ARPUS 
Мс з=  C(20Cü-2ZL Y*HRFÜS 
L. REFLECT GFF BODY POINT ON FLAME OF SYMETRY 


| 60 TOC480, 485, 460, 495, 500, 505, 510, 5152, 18 


U1KJ153=UXC1) 
Мое. 1) 
M36 1152UXC 12 
Ute ie 15zUNc 1) 
MECC Hit (SY 13 
Use dit 192 1) 
Мс +2 50260} ) 
2С 1+2 )=02 1) 
WSC 142 SUC 15 
Beco) 525 S25 431 
[S22 

1E SYMETRY 
а= усо 
bD TÜ 345 


û IFCSYH-12510, 517, 486 


MT SHER 
15-3 
2602-200 
GO TO 345 
|S=4 
YCOü--VCO 
60 TO 345 





c 


495 


496 


«Дм 


S15 


IFCSVIT-225 16,516, 406 
YZ SYMETRY 

5-2 

XCQ=-XCQ 

GO TO 343 

[S=6 

YCu=-YCO 

60 Т0 345 

IS=? 

4C0=- 200 

GO TO 345 

15-8 

REICH 


GO Ta 345 

М. ADO CONTRIBUTIONS OF ALL REFLECTIONS 
VIEH IHU TCL UN CE DUC CP FUNC CB EUS 2 
оС ЛЫ АЛ UX CS AU AUDI US) 
UCA PSUS CUP -UXCS »-UX CP 9X CE FUE CSD 
C+ PSU TCH 1t 1 UY CS SUC? SFU CE 2-1 
Е овоон) 
С + 1 VSC It AU CR 2-UV C?RUV CB 0-US CB 2 
UtcdtezozutédtTezo-Uzeieo-UzC?2MUZieMUZCR» 
2d S142 = UZE JEZ UZEE -UZ LICE UZ CSD 
USC AZUR +2 RUZ ра a UEG) FUE 
Sta PUROS 
Ua Js USE D-UXCA-UX CS 
VAC SIF PHU IC TIF PUN C4 HUY CS 
VACJ1+ 1 D=U2CJ 1+ 1 Y4 C3 
Sr AN ҢІС» 
11+2 = Ал 182)-02 042-023) 
UzCUTe-222U2C119220-U2(4 0 ТУА са, 
USC S142 SUS CI 14S 24 C4 42 03D 
ТОНЫ Оо С?» 
М2 ZO UCA 
UsCJ1252UsSCUJ124UXCZ 2 
UtcJ 1e 15s fé 14+ 15-022 
UzCat* 1ozUzod 10 122-UVC2 2 
K ы se) 
Pee ІС 142 +00762 
U2CJ 1-2 o-UZzCUI 1-225- |. 
IRIS AD У 
L=F 
UX 1CJO022UX 16. IO 2-U 16 1 DFS 1L 2 
UY 1 ouv 16 OU TCZ PTS 1CL 2 
UZ ICJ 2sUz 1C. ID 2-U TOS YES 1CL 7 
UXZ CIO 22UX2€6.100*UZ46 1 *S2€L 3 
мур = усуни со ужос > 
22 ЈО = 0226 ЈО RUZ FSA 
Wea Пот MASON USC T MRSS EL 
UY SC IO UY CAO FUSIL FSSC 7 
2361-2310 USCS IFS ) 


i CONTINUE 


М. ЕНІ ПР LOOP OLER OFF BOCY POINTS 
P=F+1 
KzK-*1 
Jz2J*20 
ШБЕК-НР 2586 
Е 4.124 1 298 
0. FREAD H 
РЕВО <04> CE 


«Е2 


K QF Е ARRAY IF NEEDED 
2412 


m m. 


— 
— 
(м 





NERBr 1593291, 290, 291 
Р. END OF COOP OVER QUADS 


999 CONTINUE 


РАСЕ = 1 
КОО ШЕШ ТО ЕЙ. 17 GO TO 525 


501 FORMATCAH FT. , ,11X, 18X, 12%, 1H, 12X,1HZ, 14X,2HUX, 11%, 2HUY, 11%, 2H2 


1, 14X,2HCP> 
2 FORMRTC?H X FLOW) 


BUS FORMATC?H Y FLUN) 


агаа еу 
a Do) mh 


€n cn Ch cn 


FORNATC?H 2 FLOW) 
5 FORMATC1H1, 1504, 10%, 15НОЕР БОПУ POINTS | 1OX,SHPAGE , 152 
PeecHix.E0.0) GO Tû 700 
МЕТЕ Ot PROE PAGE 
WRITE (5 „БОШ 7 
HR ITE (6. 6015 
LIHES 1 
LAST=33 
EC NOBP-LAST >50? 510,510 
LAST=NOBF 
DO 615 I=LINE, LAST 
FORMAT < 1%, 113, 2%, SF 13.5, 2%, 3F 13.5, 2x, F 13.5) 

û. COMPLITE PRESSURE AND EDIT 2 BASIC FLOWS 
ЕНЕ Т IRI AZ IE DE ZH ү 1 FZ 
ПОТЕЗУ СЕЗЕЛИ ХЕ РС 2PC1 A UX ICE S Uv TC UZI, CF] 
LINE=LAST+ 1 
LAST=L [| ME+54 
PAGE=SFAGE+ 1 
IFcLIHNHE-HDEF 23620 620, 
20 E v PROB, PRAE 
WHITE (6,6015 
Gu Tü GCE 
eee Sil’ 60.0) GA T 
HAITE, SIS) PROS F 
HEITE t6, 503) 
ШЕТТЕ (0 6111) 
LINE=1 
LAST=55 
ISCNOBF-LASTIT07, 710, 
LAST=NUEP 
DO 715 I=LINE,LAST 
ИЕ арта р унео) 
МЕІТЕ (5.8119 1,ХЕСІУ,УРСІЗ 2Ре РУ ЦУ ру Цео рУ UZ2C IO, CP 
LINESLAST+1 
LASTEL | HE S4 
PAGE=SFAGE+ 1 
I Fi 720 ?zü , ecu 
J WRITECG 6059 PROB, PAGE 
WRITE (6,6013 
GU TO 206 
БП ЕП 0D) GO T0 825 
WHRITECO, 0033 FAROE, FAGE 
WRITE (8,604) 
HRITE (6,801) 
LINE=1 
LAST=55 
IFCHOEP-LRST?2£07,810,210 
LAST=NOEF 
0 DO 815 I-LIHE,LRST 
E O A A ytk2HUZSC[ EFL) 
EE BA AECI ҮРСІ) 2ЕСІ) 361» ІЖәсі),У2355!12,СР3 
LINE= LAST+ 1 
LAST=L 1 NE+34 


2001 


a 
oe 
e 


114 





© 


EE 


10099 
1001 


900 


PAGESPAGE+ 1 
IF CLINE-NOBP 0820, 820, 825 


O NRITESE, E05) PROB, PAGE 


WAL TE (f, 5012 

GO TO £0 

SN 

IF «IRERD.EQ.QO» GO TO Sz? 
RERDCS,262  UX4,UN4  Uz4 

IF XEOFXS»5.NE. OQ.» GO TO 900 
pa To 920 

UX4=WS A) 

Uyg = HWS¢J+20) 

vz4 = WSCU+40) 


Гра Іңірде? 


m = = 


ЕССЕ SALL, OOO, CU 
R. COMPUTE. FOURTH FLOW AND EDIT IT 

LINES 1 

LAST=5 1 

WRITECG G05) PROB, PAGE 

FORMAT < 12H00NSET FLOW Un OTT 
FET 2S2 4HUZ z,F7.3? 

WRITE (6,531) UNA, Uy UZ4 

MENTE (6,601) 
ipencer—-LAsl Joor, S40, Si 
LAST=NÜEF 

00 243 I=LINE,LAST 

UXF --UX4*UX10L2- UY FRIA IU FU OD 
UyF'z=-UX4*UV1([ UNA FUEL] ›- en 

VFS- ZIZAK UYG EIZE | ZTE | ° 

ee UE aa T a O a 

WRITE (6,611) 1,XPCI), YECTO, ZP CIO, UXP UYP, UZF, CFS 
LINE=LAST+ 1 

BESI-LINE454 

PRGE=PAGE+ 1 

IF CLINE-HOGP 0650, SSC, Sed 

HAI TE CG 5605) FF0E  FR5E 

WRITE (6,601) 

G0 TO 855 
Ez gei 


60 TO 826 


со 


HRITEX6, 16015 JO, L, XPCJ02, VPCJOD, ZPCJQ? 
FORHHTCTBHODFF BODY POINT ,13,23H ON BOUNDARY OF QUAD ‚з; 
I SH Xz,F 12.5, 5X,2HV,F 12.5, 5X, 2H22, F 12.52 

GO Ta 530 

CONTINUE 


5. RENING TAFES AND STOP 
REMIHLD C3 
FEW THO 04 
ӘДЕ 5 
END 
FUNCTIQN SOZFOXT X2, 1, v2, 21,220 
б =. Y 4 y? 
үшүд-у 
2-21-22 
POSZERZ yE AED 
R=RESCKHRESCVY FARES CZ + 1.0E-20 
F=R+R: ¡ZE 
R=. 23FR+ ASA 
R= HtRSV/FE 
SOSF= .2Z9S*F+FS F, 
БЕТПЕН 
ENC 





RPPENDIX UI — XVZPF SECTIOH PFO6 


PROGRAM FFFOGCIMPUTZ 128, TRFE 16, OUTPUT- 128, TAFEOS, TAPEOM, 
1 TAPES= INPUT, TAPEG=O0UTPUT, TAFES=TAPEOZ, TAFE4=TAPEOS > 


Xyz POTENTIAL FLOW PROGRAM VERSION d SECTION 5 
COMPUTES VELOCITIES ANG PRESSURE COEFFICIENTS FOR 
OFF EODY STREAML INES 


OO O CO CTC) 


COMMON XPC1IOO> YPC10D), ZP¢ 1000 UX 1C TOC 0, HS C220? 

1,US 16 TOO 5, UZ 161007, 0026 100), UvzC 1002, 0224100), XEN 10019 UYS¢ 1003 

Z, VZS¢ 1003, UX8 D, UCE), VUZE), 51665 02,52(650), 536650), 41632, 0250 08 

2(3) ,PROBCISO КТК +650), 24(®5П), їх(&5б), TYCASOQ T2CG50 3 

1 (BC 1RO00 3, XT ¢ 100 3 .VTCIDO 2, zTC1002, APCS). 5144 2 ‚Sky la). 

Ir zs TUO SEXO, pix 1000, 0ч 100), OZ 100), СРС1003 TAX bru 3 

EQU VALENCE СКИ, HEC21122, CKM, МА. 2, (МР, Hz Zug CSM, HECE 

2 YS) 

INTEGER SYM |P 
С ШІ RERO INFUT 
HRITE (6,5; 
6 FORMATOS 14,94F12.357 
a FORMATOS 12.502 
sS FaRHMRTC4GHOXVZ FOTEHTIRL FLOH PFROGRRHM SECTION €, VERSION 4 > 
STE IT CIA MF 12.53 
Meee gins (41.04 StH STREAMLINES TO BE COMPUTED AT „14, 10H STEFS OF 
1 ,F3.4,22H T FOR AN ONSET VELOCITY OF ,3Ғ8.40 
21 FORMAT ¢ 1X, li i [Iis pull sr НЕТ, 5X. ЛТ ТТГ ІЛЕ); 
ЕРЕСЕПЕТС1Х,14,2Ғ12.52 
REAL COR). CFROECI2, I21, 192 
WRITE (6,909 C(FROECIS, I21, 152 

gO FORHRT THO, T2842 


C E. RERD THE FRERMETERS, T RRRAY RD SOURCE FROM ТЕРЕ 21 
READ <03) (N51), 1=1, 220) 
READS) «ТХУ, ТСТ Э, ТЕСУ, ХМС, ЧЕТУ, 2С У, ТЕСТО, е1, НРУ 
RERO (03) SKIP 
IF € HECZ20) .Е0. 2. > READEO2) SKIF 
RERO (03) (81613, 12 1, NP 
READ (03) «$2€12, I8 1, NP2 
READ <03) (3301), 1=14, MP3 
C С. READ THE E ARRAY 
HZ = NP 
HZ=éHZ+11.D 5212.0 
HE = 12 
IS = 
IF=241 


ИШ [Р = 1,} 
READ «O42 P, : 
[22124240 

li |F=1F+240 

APC i= 5 

ШЕ? - 5 

AFCA) = 1. 
НР(4з-П. 

APCS =O. 

GNC 1> = 1. 6. 

GOMA) = се 

GS? = 1.73. 

3 RERD (5,862 NOEP,NST, IEND,DT ,UX 1 , UV], U21 
риту poko HEL? [KD 
NOE=N0EF 
DO 10 1=1, МОЕ 
БЕСЕ СОЛЕ XPCIX VPCIS ZPCIHO 


CO 
r at 





(> 


i 


10 


IF <EOF<5) .EG. 0.9 GO TO 10 

NOBP=1-1 

WRITEXE,Q) NOBP, NOE 

FORMATCIHO, 15, 28H STREAMLINES SPECIFIED NOT 13) 
GÜ TO 11 

CONTINUE 

CONTINUE 

WRITEC 16) | NOEP,NST, IEND,UX I, UV 1 ,UZ I 

WRITECE, 20) NOEF, NST, DT UYI UVI UZI 

HRITECH,21) 

WEITECE, 229 «І,ХРСІЗ,ҮРСІЗ,2РСІУ», 1=1, МОЕР > 

МОЕР - МИМБЕВ ОҒ STREAMLINES TO BE TRACED. 

NST - NUMBER OF STATIONS AT WHICH STREAMLINES SHOULO BE COMPUTED. 
ПП 15 1=1, МОБР 
ХТ) = ХРА) 

VICIO = VPCID 

ZT) = 2Р(1) 

SKXC1 =O. 

SKN CI 
SK2C1) 
Te = 

IRKz5 


и н 
co 


' = 1 


P = | 

J= Í 

ОО 100 1=1, НОЕР 
P SET THE БЕТІНІ Veli ty TOU THE FREE STREAM VELOCITY 

UXicld=-1.0 

US! 16 | o0. 

MZ icn, 

UNSC | =O, 

UY2C} =- 1.0 

usu» 0.0 

O y=ü Ü 

Uz3¿4|x=- 1.0 


2 У22С1|о0-80. 


POE IPT COOR OVER THE GOURDS. 
Je. 
Fi PICK UF QUAD. INFORMATION 
X 1zBCU? 
Y IzECG Jt 15 
Ү2-В 022 
NE / EE [+2 B 
=) ر چ + ل‎ 
AECA 
E à 
XC=TXCK) 
= T(K à 
zi=TziK?} 
A =TAck? 
XX=EC I+? 2 
YREEC I+ 3 
eX-BC Ug) 
AMSEC I+ 103 
V yzBCUOUS 112 
ZVY-BCA 122 
ШЕ СПМЕПТЕ LENGTH GF SIDES OF GUAD. 
D12-502F (X1, X2, V1, V2, D. , D. ? 
BES-SUZPAXZ X3, v2, Y3,.0,.02 
D34sEQ2FOX3, X4 V3, V4, 0, 0) 
[4 125O2F X4, X 1, V4, 1, 0, 02 
ШІ CUNFUTE SLUFE BF SIDES 


R 





c 


a 
zus) 
ғат 
^ 
332 
е. я f 
uM I 


241 


“r m, 
242 


SEC 
cu 


E 39503 500, 505 
с |23-1. 

GO TO 310 
CH232CV2-V3 )/CX2-X35 
С1235=0. 

jJ = ХАТА ВН 515 
С134= 1. 

cn TO 320 
СМ34-(Ұ4-Ұ32/(Х4-Х32 
C134=0. 
Ia teed, 02 1, 320 
С141=1. 

GO TG 330 

СИ4 126v 1-4 2 /CX 17 Xd 5 
Cl4120. 


n IFCOXI-X22335,321, 335 


С112=1. 
60 TO 340 
mi =< ye 1/7 X2-2 1 > 
СІ 12-0. 
Ғ4 COMPUTE ООНПЕНРОГЕ MOMENTS 


1 ГІХАаБ +17) 


CIXV zB 182 
CIYYSEC I+ 19% 
ES COMPUTE SIN AND COS OF SLOFE ANGLE FOR EACH SIGE 
fy i2=¥2-¥1)7/012 
Me s= уз 2 723 
Сүз4сіҰү4-Ұ32/034 
DVy4 120y 174 5 /D4 1 
ЕО Тху 
CXZ3= (0 XZ-XNS N /DZE 
CX24=0X2-X4 3024 
Ci4 1=cX4-X 13704 1 
БЕ ПӨДЦЕШТЕ МА ПІНСІННІ 
EE eS Y 1, YS3,0.,0.) 
ST2=S02F(X2,X4,Y2,44,.0,.0> 
21-212 2241,242,342 
ST=312 
5. START LOG GUER THE Grr BOOY POINTS 
09 530 JO=1, NQEF 
15- 
XCOS XP IO» 
WEO- Yr CO 
AAA > 
= 
ЕЕЕ Са ЧЕ EO ZC СО 
Н. DETERMIN HETHOL 
ІЕСЕР0-57ж4 2350, 350,460 
X=IXCOH-IC PEIR HL YCO- VE IE C2CO-2C Y* 2X 
SUSAN YLO- VE EY ECZEO- ZIP ZY 
SEE AR IEC NEO VE FOR + ECO әсе НК) 
ШЕРЕП-ЕТЕ2 5 5355555, 400 
ПОЛЕТЕ ИЕП UELOCTTV. EX EXHCT BETHOU 
ласа XY VIT Zu.» 
ES SUEZ FI, AZ Y, ,Y2,2,0.2 
MS REX, 2,1, Y, 2,0.) 
BEESUPF(X X4, V V4, 2,0.2 
CLRH1zBLOGCCR 1482-0122 /CB T4P24D 12 22 
CLRZ-RLOGCCR24*RS-DzZ3 )/ CEZ*RS4U23 5) 
CLHe-RLOGGCRS4RÁ-[I24 5 / CRS«4Rá DÀ > > 
CLRA-RLOGCCRA*F 1-Dl4 1 2/ CEA4F 14D 125 
TUX-CV IZ *LLR T4 CVZSSCLRZ4CVS4*CLRSACVA CLAJ 


115 





450 


470 


gel 


TUVzCXx 12*ËLR 1+ËEX23*CLRZ+CXS4*CUR3+CX4 1EFCLR4 
TUZ=U | 

IFKABSCZ 2. 001457 3375,361,301 
сей-гжж2 

EEE SOT Yr 
Е2-250(Х-Х2 EL 
Ез-2с0(Х-Х202жж2 
Е4-2504(Х-Х4)жж2 

WIECIE TIELA 

а AL 

nee ECT X 
Н4-(Ұ-Ұ42ж55-242 
(Р«СІ1223605,202,204 


: HS 12 CCH TZ4E 1-H 1 )/CZ SR 12 


HSZ-i4CH TZ*EZ-H2 5/CZ*BE2 5 

RT IZRTRHNS 12 

ATZ=ATANCHS2 2 

TUZ=AT 1-AT2 

PEGE ZS S00, 350), 36 7 
RTS=ATANC CON2S4+E 2-H of RETR 5 5 
ATS=ATANS CCNZS*ES-HS oS C2*FiS 2D 
TUESTUZ+ATIHT4 


7 IF<CIS4 >68, SES Sta 
BATS=ZATANSCCHSIFES- AS I LER 92 


ATESATAN CCMS4 EF -He 0 CSR DD 
TVESTUZ+ATIHTO 


Reece 14102370, 370, 375 


| RT?ZRTRHCCCHA trE4-HÀ 5 / CZ FE 52 


ATS=SATANS CONE PEE 12H 12 /C 2*8 122 
TUE=TUZ+ATIHTE 
DO Tit 450 
ШЕ ЕЗИРПТЕ [HOUSED VELOCITY BY ОШНОЕНРОГЕ METHODO 
ЕРОЗ-ЕРПЖЖЗ 
BPUÜ?SCRPOZ**2 ЭЖНРГ 
Wale X?RBPUS 
XSüzX**2 
YEU 
2SU=Z r2 
FS=YSß+Z50-4, 350 
05=X50+250-4.#YS0 
HeSzzx*Gg,*PS30.*XSU»/EFO? 
ШЕз-о жүзрі ЕРІ? 
Wages AHS SAPO? 
TUAESR HS 1-C I XV AHES-C INS 2-C | YS 
He 12s /EPUS 
HdezcNecg ces ovS BPO? 
TUyzBs 1-C EXX*HSS-C]XYYIUS4- C] NEZ 
IEEE elec E CY Eee [YS OSAP OD 
ОХО = ЕХ TUSPEXVETUPZEXNCK 2 
UVCIS »-TUX VE TUAE RYTUZEVHNCR 2 
Че ТИХЕ TUE ZUERTUZTZHAKR 2 
GQ TO 470 
E ГӨНРІШЕ (HOUSED VELOCITY EY MONGFOELE HETHOLD 
НЕРОЗ-Н/(ЕрПжжа) 
WISA (XCO0O-XC 3EREFOS 
s [S= Сга тС АЕРО 
EIS (<2С0-2СОЯВЕРОЗ 
Ш СЕЕСЕЕЛ С ОЕЕТЕВВТ PRINT IH PLANE OF SYMETRY 
GO To¢4e0, 465, 400 495 S00 S08, S810, 8192, IS 
Uicuatozuxcio 
UZ =X (17 
с ИКС 





Ша АРЫ 
Y2¢J1+1 7= 1( 
USC de 1380 C12 
W1CJ1+2 =2 12 
Y2(J1+2 =2) 17 
USC S142 SZC 19 
IFCSYM 525,525,461 
(552 
XZ SYMETENW 
vC0=-Y00 
ба TO 345 
IF¢SYH- 15517,517, 485 
YY SYMETRY 


Yoou=-Voo 
GO TO 345 


495 (FCSYM-29516,516, 496 
ұр SYMETRY 
498 15-5 
хсб--хС0 
00 70 245 
500 15-5 


305 


510 


ale 


cn 
c 


іл 
r3 


vcüs-vcà 
GO To 345 
= 
200=-200 
BD To 345 
15-8 
YCO=-YCO 
GO TO 345 
M. AOD CONTRIBUTIONS OF ALL REFLECTIONS 


SEIT EU TCI NR EHRT EU CE HUN CH) 


U2GU152U26.15-UXCE ox V XS ? 3+UXC8 -UXCS 0 
DENE INES II SEUNG UCSD. 
A A EE 
2d а e DANN AUR 
MSC ІЗ ІІ НСО СОНЫ - (52 
UicJitezosutiJiTe2o-Ue2cS30-UZC? o4 u2 Ce o«Uz CS) 
Uzi Jira J2J +2 +L OS UZE? PUZ D ZUZI? 
USC 1142 DUS RUZ FUER tU CEUC 
Сот 141126425359 
M263152U2CU 15x Cd o-UX c3) 

NINE EHRE A 
Misdie122U fe IHNEN 

U2C6U0 1e 10zU2ZCU e 0RUT CA ORUM CE D 

USCJ 1e 10zUSEJ 1e 150-UV C4 o4 C3 2 

+2 1+2 09-0204 35-02637 

мас + USC I 142 2S 04 ARE 

USC S142 2=USC I 142 402 64 42 C3 D 

MICU TSU ICI TeX 2D 

E GOEN J -UXC 

VICIT ISVS EII UA ZO 

UPC d+ 19=0 10144 192) 

2 +( 1 +1 ل =( 1 + إل2ا 

U3cJ 1e 152U3601189 15-0025 

WIC A142 =U 11142 24262) 
Uz£JTez5-2U2C]1e25-U262) 

Sed 142 USC I 142 242029 

UN TCS ТОО C1 HES ICP 2 





келе" 
A 


EI 


aon 


905 


Б1 


UY 1 C UG SUYA JOAN 5 ЫР) 
Ai 
We ACI SU Z CO UZE FEZ 2 

UY2 CI SUV 26 JO 4202 »PFS2CP) 
U2zZ< N10 S22 I 4203 SZC OD 
WS CI SUS JÛ US CT ESS CFD 
UVSC IQ SUNS CU AVL MESSI D 


E JG DEUZ SCO USCS PES SCP D 

N. END OF LOOF OUER OFF BUDY POINTS 
P=P+ 1 
К=К+ 1 
Sz d 70 
[FCK-NF 2292 296,544 

Р. ENO OF LOOF OUER QUHDS 
Не-НР«ТЕКОжПТ 


00 720 | = 1,NDEF 


EEOUBHHTICZX. 13, 3F 12.5,9X.4F 12.5) 


DYE IOS- LUXLUX 0L O«UV DUX 2€ LOU? PRUX CIL 2 
DV | 3e- CUX DFUV 1C LUV DRUS 2C] XUZ PEU | 5) 
DEC) = CUXIFUZ ICL UV ERUZZC LUZ EFUZSC IL) 
IFCIRK.EQ.5? GO TO 400 

IFCIRK.EQ.4» GG TO 800 

Dü 750 [=14,NOEP 

XP EXT ODA COAH 

WE CL IEVTC 1 40VC | FH 

2РСІ?-2ТС1 926 | ЖН 
SEXCHOSSEXCIMGHCIREOFDXCIO 

SEV tI 8SEVCILOSGHCITREOFDNCIÓO 

EERE I yegEZCI M GHCIRE P D2cl) 

IRE = IRE * 1 

GQ TO 93 

H=DT 

DO 830 1-1,МОЕР 

пу CU РС е ае руно EICI D) 
D CL O8-CUX DRUSV ECL OUS DRUSZCIEOUZ UY 5) 
D2cl )e- cUX LSU 1C LUV ESUZ2CIOUZ ES U23C E) 
XFCISXTCUDSH*DXCI 3/8. «SEXCI*H 

ХТС )=ХРС1) 

UPC DOS TE ›+Н+ЖПУ С | 76. +SK | 9H 

YT SYP 

2PC (ISTO) HHADZC 1/5. +SK2¢ 1 Y*H 

Tl =F j 

SKX( | =й. 
SEY) = 


c 


5i 
Бед 
` 
ل“‎ = 
e, 
DEE 
II 
en 


IK = 5 

GO To Ge 

[FR = 1] 

DO GUS іс1,МСЕР 

OSO=OXí | REZO] PREZZO] OZ 

LF! 9=1.-050 fsa 

CONTINUE 

WEITELE, EI) ITC 

FORMATCBHO STEF, 1475 

HRITESO, 62) 

FORMAT (2x AMINE, SX, IHX TIX, TET ІШЕ Ран OX. 2HUN. 
10X, 2HUZ , 10X, 2HCE 

NECIO РСТ Ору Весту ПС EPCLID, OS 1, NOBPD 

WAITEC1G> < о р P 2Р1), =і,НОЕР 2 

ПЕТТЕ ЕО МЕТ) GO TQ 910 

ІТЕ-іТІ:% 1 

G0 TO) Sau 





ENSUSIFCIEND.EQ.O 5 BB TO 83 


REWIND 02 
REWIND 04 

ENDIFILE 16 

REWIND 16 

STOP 6 

ЕМС 

FUNCTION SOZF(X1,X2,V1,42,21,22) 
Х=Х1-Х2 

үшу1-У2 

2-2 1-22 

ROSTER ZEYH DED 

ReRES CK ERESCY HAE C2 + 1 OE-20 
Rz. 3422* (Re CRB4RE D/RO 

Re ReRE/R 

SQ2F= .2S*R+RS/R 

RETURN 


END 





APPENDIX Ull - XWZPF SECTICN PF? 





PROGRAM PFPFSTAFET, IHFUT= 125, DUTPUT= 125, TAFES= INPUT, TAPE 17, 
1TRPEGZOUTPUT , THFEOS, THFEO4 , TRPES-THFEUS, TRPE4-TRFEG4 ) 


XYZ POTENTIAL FLON PROGRAM VERSION 4 AND VERSION 5 SECTION ? 
COMPUTES VELOCITIES AND PRESSURE COEFFICIENTS FOR 
On BODY STREANLINES 


NONN 


COMMON —— X4B5585,V(558)2, ZC658), XNC550), VNCeS0? 
1, ZNC6502, UX 1£5502, US 1(6500, UZ 1C65CO, UXZCBSOO, UVZ CES) 
2, UZ2(650),UX3€6500, U/2€650), U23 (6505, XC 10658, VC TC658 D 
3, XC2 (8582, VC2CESED, XC3C6580, XCACG58 2, VCACBSED, 
4X3CB58 2, V3CBSED, Z3CBBE D, X4 C58 2, V4 (658 0, Z4 (658) 
DIMENSION XL 1500, VLCTEOD, ZLC 1500, Ux e 1580, UI C £O, Uz C150), 
ICPC 1502, GK 1€ 1502, GK2( 1502, H2C 1502, STML 150), 0885 150), МЦ804 150) 
5 OMX(650), PROBCI5), VCSC85S9, SFiS) XCRCS) YCRCS) , 
PNEP CSO WSC220) ,Х87(502,Ұ57(509,257(50) 
EQUIVALENCE CWS¢201>, NPD, (YC, C2) 
RERD(O3 XCPROBCI, I1, 18) 
RERDCDS SCHECID, 121, 2200 
RERDCOS XIX VCLY,ZCIO, XNCIO, VNCHO, ZRCI D, 
18K IP, 18 1, NP? 
REPDCOSXSKIP 
МЕБ =15(220) 
ARITECA 5) IVER 
5 FŪRMAT(4EHOXYZ POTENTIAL FLOW PROGRAM SECTION 7, VERSION | 
ІР СІМЕР.ЕП.53 READCOS> SKIP 
READS ISK TF 
REAG COD SKIP 
READ OS ISK IF 
EEROR CUXICL UV ICT) U2IC1>, l=t, NPD 
RERD(OS XUXZCIES, UY 2L 9, UZ2C E), [8 1, NP 
READS UNSC UY ZCI ERED |=), HP > 
БЕМІНІ 03 
HESCHF+ 113712 
DO SO |=1,НЕ 
ГЕМ-| ж 12 
|Б=!ЕН-11 
БЕН (042 0,(Х21523,921549,ХС26.49, 
1VC2€00, XC3CJ2, XCACJO VC4 (Jo, X365, 30, 
223642, X4 C42, 46.12, 24 CJ, CEKIP Ke 1, 2, 4818, IFH. D 
NO=0 
IFCNO.NE. | 
CONTINUE 
REWIND Jd 
DG 90 1=1, МЕ 
D 1=(CXC 101 yee2+vC 11 yee2ye1.D1 
D2-CXC2€ | женге | OA 
D3-(XC3CI жау | эжжр уж 1 Г 
D4-CXCACIH ye*z- veas | )ж®2 уж 1 С 
Q OMKCHSAMAXICET,02,02,04) 
| ҒОЕМВТ(3Е12.4,314,Е12.49 
2 FORMATCEFIZ.4, 14) 
МІПе-?5 
100 READCS, 112 UXE, WSL UZ NLIN, MAX, INRITE, MACH 
[F CEQFCS9.NE.O.> NLIN=N. 
HX JeHRXJ 
IF € HRXJ .LE. 0 LGR. MAXG .GT. NP/22 MAXJ = HF/Z2 
HINJSMID-MAXU 
MAXI=H 1 DMA 
IF (MAXI OT MIDES MAXJ=M 102 


Ж 


£0 


490 


(u 


:J GU TI 


= 


e 
na 


Ko 
(ы 





e 


E 


IF € HINJ .LT. 19. MINJ = 1 

HRITEC? ?  NLIN 

HMRITEC172  NLIN,UXI UIT, UI 

IFONLIN.LE.0) GO TO 550 

МБІТЕ(6,30) «РВОВС12,11,159 

HRITECB, 34» UXI,UVI ,UZ I NEIN, MXJ, [WRITE AMACH 

34 FORMATCS4HO0N BODY STREANLINES - INFUT DATA ¿6H VXI =,F 

16H UV| z,F10.5/6H UZI z,F10.5/5H NLINs, ! 107 
2 6H JNAX=,110,/,SH IMEITEZ, 110, 7, 6H MRCH Nü-,F 10.5» 
HRITECB, 38) 


Se FÜRNATSZTHOSTRERML ING STARTING PÜIHTSSA LINE, 11%, HA, 123 


ШЕГІ ІНЕ 10X, S<HNSP A 
LIN=NLIN 
[ü 45 I=, LIN 
I len 2з Х|» БТ), ТС, МЕРІ 
CE EOFS EG. O.D СО Та 45 
NLIN=I-1 
HRITEGB, 42». NLIM,LIH 
EE SAAT: ROU, TS, ZEA STREAMLINES SPECIFIED HOT ‚| 
ЕЧЕН. СЕ. С> 600 Тп 550 
GO TO 48 
45 HRITEXO,925 1,8 
4A FORMATCIX, T3 Zx, 
4e CONTINUE 
Мей PEREZ | ee? FES 
ШИ РИНЕН Еа И ç GO TG 1130 
ECLOSUBHFUTE CRITICRL MACH NG. 
USD = OQ. 
DG 1100 J=1,NF 
US = ARTE IH PRIZE Ue DEUS Cf Sebo + 


“ә 


б) 
> 


* тое s = i б. ` 
Dee ea el el 


TOL}, 
F13.5,19) 


لر 


í. 


1 ОЕЕО ЕСЕ УОТ урен + 
2 tous Toy Nee р тел мел р лк 


IF (65 Бї. USO) USO = US 
1100 DCOMTIMUE 
U = SORTCUSO SUSO > 
man = 1/1 
Ш 110 isi, 
Mi = СОСЕМА ои yp Pos 
СИН Чек сы ЗБ эж] ТУД 
ШЕТІН zs CHHRSCHHC-CMNHESSe? Y/CCHHRECMHE-Z SECMBE 5 
KRITELS, 11203 CHHR 
ШЕН ЕПБИНІСОІН CRITICAL MACH H3. z,F5.3 
moa CONTINUE 
START LOÜF QUER STRERHLINES 
me 300 =], М 
CIRT=1. 
161 dJi=1 
AF=1. 
ШУМЕ =й. 
Jy ceM DOD. 
UECHIO =0. 
СРЕНІНШПО-О, 
HMDS t. 
GE TCM oan. 
GK24MI ED 220. 
STHL CH IU =O. 
102 Nü-HgPCLL?» 
LHüsHU 
BECHIBOZXSTELL > 
VLGHIDOZVSTCLL 5 
ZLIMIDI=2ST(LL)D 
dati 


7 


124 


10.94 


t, 1HV, 


“sy 





cy 09 


C) 


105 


600 


610 


ball 


107 


ШЕЛ 

SEPARATE CALCULATION OF SECOND 

POINT FROM MAIN LOF 
XLT=SCXLEII-XENDODIEXSCND IH CN LEI I-VENGIIRYSCNO ) 
1 Erz CHI eZ SCHED 
VLTZCXLGUIO-XANU 2 PFXACHQUHE CYL CU 0- Y CHO J EYAR 
1 + (ELIZONDO S424 CHO D 
XL DAL TIA CNO FY L TEX CHO SAN CNO 2 

VLC J EXLTPVYS ONG D4VLT#¥4 ENO DEY CNG 9 

ZL CI EXLT#Z3 (NO IH YLT*Z4cNO ^*-z (HU ? 
IüTZMODCAHO,42 4 1 

GO TO ces0,600,610,6205  IQT 

NR=NC+ 1 

I= 

GO TO 107 

NA=NO+2 

NU=NO- 1 

GO TO 10? 


1 МНаМП-2 


МЦ={Чй+ 1 

бо To 107 

NFizHU- 1 

NU=NÜ-2 

A IU ECHO DEU ІХ Ма 92 

UYO=- CU PRUY 1 NG OA DI UV ZHOU RUS CNO AO 

Li2G=— CUNT AZ PONE TEUZ CNU | AUS CN 2D 

ЕЕ TCHR UY AUX CNR te PEUX SCONE 2D 

URS- CUN EUY EHR AU EUY ZLAR OUZ | HUY SCONE 3D 

UZR== CUP RUZ CNR DTU PUZZ CNR AUZ | HUSS CME DD 

IXUS- CUX [UPC TCHNLE HEUS DRUUCPCHU OU РЕСМИ > 

UN Liz - CUN peu? T ONU SUP DEUSZ CNUDHEUZ EUY KHI) 
Eus-cCD DRUZTCHUO UM pRUZZCHU X Uz LELIZSCHUP2 
TRANSFORM VELOCITIES TO üUHL SYSTEM 
UDzUXOTXSCNO MLINVEPEVS CHO 26 FC Z SCHO » 
VOUT ONO FURY CHOOFUZÜFZ4CNO 2 
CSRES1./CXBHCNO XEXNCNBE УМЕН ЖУНСМЕ 5+2МЕНО ЗЕМЕМЕ 5 2 
UT=SUXFHAS CHF IH VEN SCHFR IH IZRFZICHR > 

UT= CURE CMF LIM Feb CN FU ZA 4 CNL 3 ЕСЕ 

XXF= CXINA EXCH НИН ОЖ На ЭҒ2 ЗІНЕ 232 
AYR= CX4 CHA EXONS 244 CNA VS CND 0424 CHE YES CNT DD 
URZUT*XXR-TUT*XVR 

UXHT CXS CNR PEE CHD OVS CNR FYE CHEF 2S CHR MSE CND DD 
YVR= CE CHA PERS CHO t's ON OEY SOG 423 CNR OSS CNG DD 
URSUT EY XR SU TAY 
Uu=UXUFAS CNO AUS CNO tU FZS и, 

MIEL OXON CHO ACH CHL OF YN CNC жун ОНЫН EEN CLL Oo 
USC UU CM SELLY CHS FLIZ Lid CN D8 SOSH 

FINO RELATIVE COORDINATES OF NEIGHESRING QUADS 
XD=X¿<HF I-XCND > 

YD= (MAA У 

202 CHA I-2 CHO) 

XT SXO#KS (NF +N O'S CONF 420423 (NF > 

YTT2xU* x ANE ot Der CF 2+ 2024 HF 2 

21 XDTXHONP 2E VDEENHE CHE 2e ZDF ZN CHIP, 2 


VTeC-ANSORTC/TTES242 T9x2 )2  TT*CSRe VTT CER , T6666? 


XF=XT*XXF+VT*XVF 
MAA HEYF, 

AD=XCNLII-X CHD D 

YESSY CNU J-Y CHEO 

eD-22zinu -ZNO > 

AUS LOS ONO ER 3 ONO + 2D 23 CONO) 
YT= 1044 CNO y+V[PEV4 (HG 3+2Z[rF24 (HQ 2 


25 





putt 


ZTEXD*XNCNO Ye VOPFVNCNO 20 *ZHCNO ) 
Us (4 "SORT CUTPE242T9*2 VT /CSLI& VT »* , 16556587 
FING! COEFFICIENTS OF VELOCITY FUNCTIONS 
DEN=1. /(XPEVU-XLPYR 9 

Uu 12 CCUE-UD YRLI- CUL UC Y*V REO DEN 
U2--ccug-Ua*Xu-cuü-Ua XR »*DEN 

Ute CCOURI-UO LI CULO DAYA) DEN 
U2z-((UR-UD XEXU- CUU-UQ »*x&. EDEN 

FIND VELOCITY AT STREAMLINE POINT 

USL=UD+U 1XLT+U2*YLT 

USL-UOU FEXLTUZSVLT 

UXPeLiS X3 C NO +USLAYA CNO) 
LIVPSLISLAVSCNO USL CHD 
UPPSUSLAZICNO USL #24 CNC 

FIND GEODESIC CURVATURES GK1,GK2 
VSOD=LISLERZHUSL HH 

DENSUSOU*SORTCUSOD) 

DK IPS CLIEL e CUSL*U2-LISLFU2 USA USA -USLE IDEN 
GHOP= CUSLECUISL AL USL'*UZ -USL#CUSL EY USE OO ¿DEN 
FINO LOCAL STREAM FUNCTION 

CXS CU HUCI- LIIE USCC 
сүч=Ш2-ий 5 ДЕПП 

СХХ=Ш2-СҮҮ-\1 
COSXLT*UO-VLT*UO-CXVEXLTSSVLT-CUVevL Tes2- CXXXL Tz 
FINO STREAM FUNCTION AT CORNER POINTS 
XCFC1=xXC CHO 

XERKZISKCHENG) 

XCRCS2XC3CNO) 

XCRCL oEXCACNDO 

MCACS ENCORE 1? 

WAY D=¥C1CHOD 

WORC2 =VC2CNOD 

VCECS 28 VCSCNO 

VCRCA ve VEA CHO 

WCRCS ISYCRC 1) 

00 110 N=1,4 

SF LOS COYOACAN IFLR YCR CN DFO CECH ICR CH EVY VCR CN ож 
I HEXXEXCREN JEEZ 

SF¿S)3=SF(1) 

ШЕСТ-0 

DO 120 N=1,4 

IF ¢ SFENMSFEN+TI GE. Q. > GO TQ 120 
XMEOGCRODAXCRONA 122,5 
VMECUCRON S VCRCNE 12.5 

FIND) INTERSECTION WITH SIDE OF QUAL. 

SE MECC IAEA CER pp OEY [kg 
RC=2 *(SF¿N3-2. *SFM+SFZH+1)5) 
EC=3.*SFCN)-4.3SFMHSFCH+ 1) 

[Е сас £9. 0) GO TO 112 

SRESORTC(EBD**2-4 FADASECNY ) 
TF-GCAHSRO/C2 HAC 

IF «TF .LE. 1. .RMD. TP .GE. D. 2 GO TO 115 
TPz(EC-ERO/CR AC) : 

60 TO 115 

IF «BC .EQ. Q0 » Gü TO 120 

ТР=Е (НЕС 

ҮНЕ-(1.-ТРОЖУСВОМОТРЖУСВОМЫ 12 
VAP=C1.-TP YEVYCRCN + TRAYOR ONS 1) 

TESTP=¢ CXNP=XLT DFC QVHP-VLT »*UQ ED TAT 

IF CTESTF .LE. TEST) GOTO 120 

TEST=TESTP 

XNT=XNF 





YNT=YNF 
6 CONTINUE 
IFC TEST .E4. 42 GO TO 2€0 

С AVERAGE LAST VELOCITY AND CURVATURE 
l ОХИР SEAS 
UV CUI = CUY Ed ALIVE ERE 
UZ dos QUZiJo«UzP ARA 
GK 1C J5oz CGK 16004 BR TP. Y* RF 
GK2¢ Jo=CGKSC J AGREE »*RF 
H2C Jo=HoCUL o#C 2. -GK ICUL HEC STML CDO-STML CL Oo oft? tb RT 
1 SMOC 
CR I. - RC ERE ZERO жаа FZ DNS 
MES OIE SAA TOY. CRD) 

C Bun UTE UEPDCITY RT NEXT FÜINT 
NOLA AHA 
«Дт. 
J=J+J| 
USL= IS ANTALYA TL 
ESESUOTXNTRUJ VT ez 
UC J SUS LAS NO IHUSL EX4 CND? 
L<. SUELA CORO MUSA RE 
Пё SUELA CON A SAO 

С COMPUTE GEODESIC CURURTLAGS 
USODSLISL*#DHUSL Y2 
DEHZUZD*ZCRTCUSOD 2 


a 
to 
С. 


бК2 De CUSERCUELRU He USLHU )- жп! жн жи2 Уу/бЕН 
COROSEORTCCXNT-XLT ORZ4CYNT-WLT FZ ECL FIT 
ЕТМ СТМ С + CORD 
C —— COMPUTE H2 
Hz. -HZ CL YeC2  -CORDGE 10.L 2/02 &CORCYGE TX DD 
CPC 1, -UEen /U8n 
URESCJ)-SURTCT.-CE CJ? 
AF=.5 
LH0=NQ 
XL CJ =XNT*X2(CNO D+VNT*X4 (CN FY CNC 
LOSA TAS CNO MH VT COH CHO) 
ZL GIXSXNT*ZS QNO V HTFZ4 CNODZ CHO 2 
IF € J ДЕ. MINJ бв. А GE. NAM 3? — GO TO 280 
C FINO NEXT QUAD. 
le 
250 NOI 
IFCL.EG.LNG) GO TO 250 
TESTO Y RL E] ye 
1(2LGI-zC IE PPRRZ-OHZEN > 
IFCTEST.GT.O.) GO TO 280 
A FZ 
A ее? 
DSS=CXCSC | INCH S| EEE CVES CL )- VC4 CI EE 
d= (NCIC LINC ACL MRDECHLA CL HE TEL #2 
MLTROXL EG IEXNCL EXEC CYL CIAL жүзі 
Тар арусе 
YT OLI ACA GD- VCI Y Pd C | >+ 
MELO ZC ZOO 
ZLTe XL GJD-X€ | 338XNCI E CL GD-VCI YHE I + 
MELO ZU IAEA) 
250=2LT++2 
TEST=250-. PDMXC1) 
IFCTEST.GT.O.) GO TO 280 
RC I-SDRTCZEQe YLT-XC 1€ L2 Ye*24 CLT-VC ACE) 
pez ) 
RC2-SORTCZSDe OIL T-XCC 1 2 Ye24 CVLT-VCZC I 22992) 


1 


(0 7T 


< 


127 





ЕСЗ-ОПВЕТ (2604 (21Т-Х035 оя ONL T-VCG ESD 
RC4H=SORTCZS0+ CXL THAC SHC | oD EERE CYL T-Y4 Cf RFD 
TEST= (¿RC I+FC2 #2 U-019 1 *4 21 
[ЕСТЕ Т.ЬТ.П.) бй Тї 105 
E TEC TELS y=DSZE *1 21 
IFCTEST.LT. 0. 1 Бо то 105 
TEST= ««НВСЗФВС4оОжж2;-П63 31.21 
MES SET.O.5 BD 10 10% 
TEST= cCRCA«RC 1 »**75-D54 *1.21 
J TES IIT у со тоа 105 
250 |=1+1 
ТЕСІЛЕ NP y 50 TO 250 
ЕШБІР ПЕТ LT. OG.» GO TO 28 
Піріс- 1. 
Jl =- 1 
JAM. 
60 TO 102 
235 JUMIN=J 
SS=STHMLCMIND 
DQ 290 J=JHIN, JHAX 
СШ ЕТМЕ < C J5=STHL ¿J y-SS 
J= N+ i 
JAS IMAZ 
r=}. 
к= ч 
ЕРЛЕР ЭПӘКРРШЕСІ>, ісі, 15) 
gl FORMATE THT, 11445 
IIE, EUSA UYT UZI 
IEE ТЕЙ CREE EEGOH. UXIZ FO. S,2X,4HUYIT,F6.3, 2X, dHUZIs, F65.32 
ВЕШ ЕЕ оа DE МӘРЕГЕ), XSTCLL 23 TULL сы» 
SQ FORMATS THO LINE Nit. pee S PASS ING THROUGH QUADRILATERAL  ,13, 
26H WITH STARTING POINT, E STI IZ4 E. 
2 A 27) 
EN UUHCEE A MINI HE. мА BE MALO ES 
ERANS t SoH FPROBARELE ERROR = LINE 15 VERY LONG) 
i ER 
fee, ES STML¢J+2 -STMLCL- ШШЕ ЕС жге ЕСТ ELD) GUO TÜ 320 
КОШЕ СІП), XECE OU VECL гісі), ir EET ПЕТ СІ; 
ШЕ: КАП К ӘП ЕШІНІ DELETED ,10Х,-Ғ 12. E Ed ch 
STMLC OSCBESSTHL COL OSTHLGI DOCE. 5 
ХісісзсЕНҢЕЖУІІ Сі ох Сы ОЗ ЖЕЕЗІ. 2 
AAA A dd O A A t 1. 
SL ARE EE LEA IA НЕ]. 
IIX<L3=SHF HL CL IH LI SH I I I / CHE * 1. 
MESS A CL dr Nr III FIAF+ 1.) 
LIEILI=SCHFFLEIL HUF IHN EREI. 
GE ICL G=CRF GK PCL A4+GkK 104 1 oa SCR +1. i 
biz. OSCRETDR2CL HF GER ZGUM 1597 CRF 1, 
На ISCHFFHZIL IH HF d12027CHF4 1.2 


^ 
2n 


Si 


оз 
— 
v 


Г. 5 UA 


ЕЕ ОХЕ очу СЕ ужин С oe? » USO 
MES EL O=SOA TOA. =CPEL IO) 

AF=AF + 1. 

GO TO 520 


' ТА 
ко 
c 


RF=1. 

L=L+1 

Kalt | 
STHL 4L 2-2STHML CE 
XL CL )eXL CK 
YLCL Om VL CKÓ 

eL cL »-2L CK 5 
UXEL ISUXCKD 

UY CL =UV (K) 
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31 


Ë ЖАС 


1141 


1150 
1150 


o3 


50 


e 


Educ 

GK 1CE 28 6K 1CK 2 

GK2CL >=GK2<çK 

HACL 2=H2 CK 2 

CFL 2=CP(k 2 
ABS CL =URESCK 5 
NOLIRO £L 32 NOLIBRD EK 2 
COHTINUE 

L=L +1 

STML CL y=STHL < JIHRX 2 

ЕСЕ = СМЯХ > 

VL CL 22 VL C RX 2 
ісігі МН 2 

Ux CL 22UX € JHHRX 2 

ШУ =U CIHAN 2 

ОЕ = ЛАХ?» 

GK ICL =GK TC UNAM D 

BK2 CL d=GK2 ¢ IMAX D 

Hz&L 22H24 .IMRX ? 
CELIS EF CUNA > 

URES L SURES 4 IMAX 2 

JIMASEL 

NORAD UNAX 3S HOLIBUC rax - 12 
MOURSDCUMIEN SNOLIBSDE UM HR 12 
HRITECB, S512 

FORMAT C4HOG — 1, 5X, 1H, 9X, 1HY, 


19%, 182, 09%, 2H, 8X, 2HUW" 8X, 2HUZ,09x, 
22HCP, 8X,2HK1, 8X,2HK2, 8X 2HH2,8X, 2H5L,8X, 1HU,QX, 1HF) 


Heine ЕС GO.» аот 1160 
COMPLITE COMPRESSIENLITY CORRECTION 
D0 1150 J=JMIN, JMAX 
CADA AD 
USDA = USD 
SM = АМАСНЖЖ2 
ПО 1140 1-1,3 
Ff = (4.4. 2#5M#C1.-USDR DD 
ECE LT. 000001) B = 000001 
VSGE = USD ¿ARSS 
R = 1.4. 29501. -NSTE DD 
IF CR LT. 0000012 Б = .000001 
- USD/R**2.5 
VEDA = CUSDOFUSDA-USIR#*2 > /CUSDC+USOR-2 . #USOE > 
C14 DESM#C 1 -USOAD) 
IF CR LT. .000001> Е = 000001 
Б = Б=®] 25 
ШОО! Уу = ИСЛ 
ШОУ = WEHR 
Uz(J5 » UECJD/R 
URES(J) - SORTCIUSDR? 
EP£UdO 2 CR**2, 8-1, 7C. 2*8HD 
CONT 1 HUE 
38 
SS 1=JMIN, MAX 
K=K+1 
HRITECG6,505 K,XLCID, VECID, ZLCI XET UFE UZETO, CPCI O, 
GK1CIO, GK2CIO, H2C IO, STHL CIO, URBSCIO, Natenic po 
FORMRTC1X, 13, SF 10.5, 1X, SF 10.5, 1X,6F 10.5, 15? | 
FORMATOZF 12.5) 
HEITEZ172. K, OLCD,VLCID,ZLCID, NOUBDE ID, ImdMIN, JMRX 2 
IHRITE .LE. ü -- WRITE SL,U,HZ,K2 
IHRITE .GE. 2 -- HRITE X, V, 2, CF 
IHEITE .Eü. 1 -- WRITE SL,U,H2,K2 RND X, V, 2, CF 


"n n u 





IF CINRITE.GT. 1) GO TO 340 
WRITEC?> K,  CSTMLCHD, URBSCH), HZUIDGK2CH), T=JMIN, MAX) 
IF CIHRITE.LT. 12 Gü TO 400 

HRITEC?) — K, CXLCDD, ECDD, 2LCIO, CP£IO,. I=JMIN, MAX) 

GO TO 400 

NEITECO,S0) NSPELLO 

NRITEXE, 65) 

GO TO 282 

CONTINUE 

Go TO 100 

READ NEXT SET CF STREAMLINES 

ЫРІТЕ(6,451216,М0 

FORMAT ¢ 14H TAFE 04 EFFOF 214) 

ENOFILE 7 

REWIND 7 

ENDFILE 17 

REWIND 17 

REHIND 04 

STOF 7 

ENC 





RPPEHDIX Uill 





SAMPLE PRUELEM TRIAXIAL ELLIPSOID 


280 


1 


S 150 


T - "- m 


оосо 
97861 
‚90789 
asco 
‚21563 
,62478 


estas 


Ms! 


‚44721 


AIS р 


Ns sul 


Es 
. 23532 
. 15957 
. 11457 
15245 
КЕШЕНІ 
ا‎ 


TT 


Pac 


GEE 


=" 
c 1 ' 


‚Сыды! 
221494 
ES 


. 9247 1 


Sq Fe 
бата OT 
RES 


т лы 
Leo 


ж m = = = 


fils 


. de'm = ` 


. 16346 


B 


‚й524 1 


‚99499 


¬< 
d lof 


00234 
$1347 


“яла 


zent 


£ ^40 4 


QUE TUM 


dez dd 


y ee 


44497 


po 45 


> 
еер ce 


20672 


әле 


N 
16552 


. 11410 


neon 


1 
MERCI. at at 


‚00000 


” 


Fur a4 
. 89702 
281422 
0004 
.6 1771 
521843 
44215 
EST it 
.2u48 

‚22424 


Eur 


50 3 


ПОГІПГІ 
, 000E 
.Qnann 
00000 
00000 
. 00000 
. Qno 
. 00000 
20000 
.00000 
. 00000 
HONG 
‚ad 
. GOOD 
топо 
. Poe 
. 10000 
. 10900 
. 10000 
. 10000 
AGGI 
. 10000 
ШЕНІНЕ 


ДП 


ma = ‘= 


12000 


. Taco 


. 10000 
10900 


AE 


P 


20909 
. 20090 
220000 
¿20000 
. 2000 
. 20660 
205022 
20000 
2008 
CUI 
. LOGIC 
20000 
. 20000 
2010010 


‚эпте 
220000 
BCEENEK 
20000 
30000 
¿30000 
UIO 
¿SOCIO 
. 300100 


00001 ü 
. GOCÖC 
. 10250 
. 20980 


: 22359 
. 340 14 
390140 
‚42544 
‚4472 

46243 
.47Ү25 
42276 
22275 
40570 
49331 
ОСС 
00000 
‚10275 


20924 


ғыз em mes, 


5 айласы ҒАС, 
24270 
. 38991 
4240 
44555 
‚3257 
47855 
‚42515 
492 15 
2490808 


49359 
‚34437? 


ССС 
‚10254 


arose 
ella 


30347 


‚ы 
34739 


тост 
. œh in 


ut 


42330 
„44437 
35113 
47425 
. 
45029 
.4942 1 
4985 4 
498749 
O00 
10184 


- TRIAXIAL ELLIPSOID 


+ Pp. P. fo b Ahh OOOO WH) оу бз о С) С RIND RD For F RO NO Ba RO FO RO NO FO RO ca — — — Là La — — — — — — — — u (СӘ 


INPUT FILE 


0 0 -3 HN LA 60 H3 — 60D 


ہہ ال٣‏ ره 


C 


mh u li «т=з ==, 


ua CO -J Ch £n La 69 ho — Cn R09 hà — Deck CO =) сл сп Pa 


m" 
-— ےا‎ 


2 


к 


” 
= 


cn E 


— — mh 
F — C O CO —J en CN La 65 S — 


— — — — — — — — — — — — — — — + — — „д — — هت‎ — — — — — .— — — — — — — c — — — — — — — — — — - „д — e — e o o — í 


ity rJ 


С Сз су гэ сс а Су С ЕИЗ 


EN GC SS 


"OD O ere 


CN NET ENN IIS 


CODO UD OCC COTO O 


O 


Ë 


SEI 


‚байт 
СОСО 
 acono 
‚00990 
. ODODO 
. Qüpnnü 


Good 
00000 
‚аабай 
0010100 
00000 
Ba 


. 00009 


= «m = © == 


Qna 


CORTO 
000110 


СИИИ 


СОС 
00000 


ЕЛЕНЕ 


. 00209 
‚059000 


nn m 


СИС 
001110 


‚босо 


= - - d "A 


GG 
(00220 
JPI 


. 


mm “ав - “- 


OOOO 
‚болоо 
.000090 
(Шой 
GOO 
.O0000 
Qon 


‚Ой 


‚Заб 


Xn 





(алан Сә -а (2 


C 
f 


n CG CG C c со 
1. С2 Саар Co 


n 


ою шоо о n 
05 (I -J € Él — 


‚йезе 
ШІНШІ? 
EES 
‚51476 
43818 
28771 
РЕСЕ 
23214 
‚ 16624 


11235 


‚05142 


Сино) 

Grua 
‚9522.4 
SEUS 


270137 
ata 
‚51476 
‚32312 
ESE 

‚23219 
22214 
. 15624 
ШЕР, 
Coa? 
. 00000 
mus» 
84754 


Seeds 


mee tt 
OSCE 1 
QOO 
‚95334 
43354 
. SED 
. ESAD 
ЕШТЕН 
‚93500 
seit? 
42562 
ESS5U 1 
20448 
. 22601 
mG tao 


10820 


. 30000 
. 20000 
¿30000 
OCIO 
.40000 
. 40000 
. 40000 
. 40000 
. 40000 
. 40000 
. 40000 
40000 


- m = 


. AQ ul 
40000 
40000 
40000 
40000 
. 40000 
40000 
¿40009 
40000 
. 4n 


F000 
40022 
40000 


ADDIG 
, 40000 
. 400029 
OOO 
400d 


ИЕН 


. 


gana 


40000 
ЕЛІГІП 


ccd 
. SüLCIOD 


. aci 


50000 


STA 


соо 


. SQ000 


СОС 


ЕЛЕНЕ 


done 
. 90009 


.БОППО 
. Босоо 
enano 


‚бойоо 


ЖЕНЕ 


60000 
60009 
‚60000 


ADODO 
ecco 
¿GOO 
Ма ИЩ 


.49 108 
IDO 


. 49454 
кана 
. dac? 8 
eium 
27780 
.a4208 
‚35251 


. 4 1554 


.43818 


. 45409 
.467ғ1 


LIN. 


48280 
4888? 
48922 


. 42990 


OOOO 


¿OS 


aa 
. 


TEE 
54202 


ANNE 


41654 
‚423318 
45409 
48761 


A = 
: 4 vo шо 


48280 
3566607 


дє" 
. 


. 


CIDO 


00459 


20295 


+ 


La a fa La Eb Fa Jas (5 


tt am em 


сті 


ші) 145 (Га) 1. a Sint A 


=J O لا‎ се Q) -a DA cu 
tow G3) CO -J CO 1) 


Fr —J o —J Ch Pa Qt a ol) г. 


42 


48346 
48412 
‚00000 
‚033 12 


en S a 


ес 
239242 
.+0554 
42681 
244211 
49920 
48250 


70105 
тасу”? 
ss. 


Cn Ch Ch Ch Cn 0i Cn Cn cn 7c cen Ch CA CA CA eon «oun en en en en ¿n CA CN en en cen cn p A p+ 


ст 


QO CQ aj a) a) =J =) ل ل-‎ —J —1 —J —J ل“ فل- ل“‎ Ad آ)]‎ 


^ CO co ue 


го түз Q3 CO сои со чуу гуп 


mm y (REN AI — Cn -R. (a 


a O iD on -J Cn ¿n D. có MN — 


— 


— 
23 hà 


i£ Qo -J Ch CD p. 6) à — CH P 


4J Ch CA fe LA pA ع‎ 


КУ МУМУ МУ КУ РУ КУ FE Po fO FO SO FO fJ RO FA FO RÛ SO КУРӘ RA SO P RO FO NO HO PO RO Fo PI BJ MK RO BO PF) RB FP P NO FQ — — — — — — — — — — — — — — — — — — 


сл Сз гл Сг су Сус Са СМС СЭС Cy ODG AE 


— 0 


Са Сес (ка = Cw v eS 


c 


OC E O 


MIA) ES ша ша кн =i 


С» 


C cC) 


800200 


гай! 
00000 
„сабаб 


. LUNI 


0012100 
00000 
. 00000 
00000 


. 


. 00000 
¿00000 


= =, =, m 


КЕНЕ 


ж A am ә 


асс 


‚00009 


000 


‚Good 
. (O00391 
Ообо 


= ~ = 


AGAL 
. 0000n 


ODDO 


ПІПГІГІГ 


-G0 


ЕДІН 


QADA 
Сес 


UCU 


=, ت‎ =, а=, 


SOU 


Телеу жеде, 


ею - m бев 


‚22000 
‚00000 


оса 
жасай 


CARO 


Coon 


e "`` % s m 


ET 
.00000 


а» е = o аә 


Сеп 


‚Ойбай 
Con 


eo 
ШЕНІНДЕ 





05006 


=.” 4 rr. 


a E ta 
EC Sc gc 


. 60000 


en. 
‚Lilli 


4°53 1 


. 192 10 


11 


QUOC 


ncc 


E i4 2 Û 

‚Dich 60000 476097 8 15 2 A . 90000 
. . O00 ‚ooooG 4 1 Ze асап 
‚0161 70000 а шр. 2) 2) O 0000 
¿3047 FOCO тоол ал Ü Onion 
27150 ‚70000 тере dard 2 Ü nani 
¿67056 20000 [сс 757270 ‚Обаа 
Paseo . OOOO Rosen H5 2 Q 00008 
.492 14 70000 ssc ea т > Ü (00000 
418423 PORO 41993 9 8 2 Dn ЖЕБЕ 
. 3a 156 PACO 42414 а 8 2 0 anna 
227436 . 0000 Salem a п 2б Иза 
222144 FTE 45504 9 11 г . aco 
15894 700100 ISE N 00000 
10742 тоссо 46579 9 13 2 б -00000 
04216 . 70000 36772 8 14 2 (0 80200 
00000 . 20000 46537 4 15 2 Û 00000 
21832 ‚sooo (50000010 o d 00000 
804801 кана C9477 10 2 Dn , OOOO 
23210 . SONO ‚19210 10 24 £00000 
275424 ¿Sauna ESI. 10 Ze , OOO 
65507 EDO 31999 10 Ze abl sua riti 
97262 .80000 235781 10 PESE Goch 
„46151 ‚Sag ros UE то оч Oriana 
40952 80000 40858 10 ШЕП Ono 
34347 . ador 42478 10 DE ПІРІН 
237 solo TAP iO 1 Z OC 
21714 zn 24452 tu 1 A 0 борап 

Saal eut 41810001 2 Č Croce 
iz Lu Ca Chai чол TE 1 SE nonni 
04310 есеюі 45763 10 1 2 .Bccao 
GOOG OOD 45226 10 ? 2 00990 
a1557 25000 ¿OOO 11 2 ¿00050 
‚2989 1 30000 ‚00427 11 Э Өй 

a 
3 


ы олһз-- ASA SS SES ED SN AS nn CH bh O N iN ا‎ whi Ra - A p су М — C 


275454 Silay 25481 11 aa 
658507? Scoach „21999 11 OCOD 
Gro? eodd „зата 11 unco 
.38151 Сн ‚seas2 11 ОСС 
anas ead 40988 11 e Coen 
34297 80000 42478 11 3 ‚00000 
27322 ennnn 42741 11 1 = aoo 
21714 ZATEN (445210111 = ‚coat 
15950 .Sanaü .48151 1i 1 3 Dua 
METIO 60000 ааа ТТТ 3 З . Qno 
04810 ESA 5768S dii 1 3 Кана 
к ШҮ SOC S226 11 1 3 00000 
¿29303 «абас! Nu SEI 3 . 09000 
‚87333 02000 ‚09155 12 C5 .00000 
КЛЕР? SOCIO 18718 12 3 Ü HH 
7590 . SG Bar Ze s Ü CUO 
52926 S00 Sfl 12 з д ‚бип 
285744 . aco 24554 12 SN D OOH 
459017 30000 209959 12 c i ¿OOO 
. 39037 anodoao ¿39937 12 з 0 90000 
32515 asada 41388 12 S 021000 
5632 а0000 2520 12 1 з Û . GOGOO 
ЕЗІ ЕЕ anonn A2290 T2 1 з D pagni 
. 15152 . UU 44004 12 1 са о ‚ООо 
. 1024 1 . aaücn 44357 12 1 з 0 . aac 
¿04637 .90000 44590 12 1 so Û S 





‚00000 
. 66603 
84780 
. 48625 
. I320 
.6 1995 
2410? 
‚45495 
387320 
ала 2 
асар 
. 208 182 
‚ 14®е4 
09931 
04545 
О пор 
. 33316 
21730 
¿20329 
¿Brad 
29794 
E 


m 
43577 


Ese 
[51545 
.24u0n5 
. 19727 
. 14179 
UST? 
‚04283 
. Boon 


econn 


-— tu 


NG 


— — P каз a Тал CH ch 
“зго G3 (D cmn r.a uy 4 LON Pa 00 


' 
С 


sJ -J 90 CN ao OD Cn 605 £o 00 


0ا س CO‏ لہ ل۴٣‏ ت = حا له عط لا A‏ 


tun o c) 8J Cu SO CO vn pad 


бай 


шаи 
. ec 


IO PAG Cf pa 
A (бз CD 


ı BJ -JI CG — 1 


RB cn c 1 oJ 
243 -J ON pad) 


end بے کے‎ omuia oousi h od h os h h k h als h h os rm e ls شس‎ dl los h شس شس‎ s h شس شس‎ f h cnn 


LU E P ть شس‎ 


بش rs‏ شس شب شے شس „л cenae‏ شس he м s...‏ شس — 


— 


LH E M P ee eth, 


(80000 
ООС 


. 00000 


. 00000 
00000 


. 00000 
00000 
.0С000 
. 00000 
00000 
cag 
. 00000 


OCOT 


. пайса 


20050 


on nm 


. 10000 
. 10000 
. ONO 


. OOD 
Woo 
. 1000 
. 100800 
. 10000 


. 10000 


. 10000 
. 10000 
‚отт 
TOGO 
. 10000 
¿AUTO 
20000 


SOL HH I 


. 20000 
. 20060 


2001010 


MOD 


о de = œ œw 


. 20005 
. 20000 
РІНЕН 


. 20000 


. 201006 


РІНШЕШЕ 


. 2005 


20000 
20000 
20000 
20000 


20000 


. 20000 


.2ü000 


. 20000 


20000 
. 20000 
. 20000 
. 20000 


“2 
. Kraul 


^A 
И 

теле, 
221 


бал O AA 


бу Су Су Сл‏ ج 
(uie b‏ 


. 44651 
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